Goal: To determine how chronic infection alters the ability to survive secondary infections of unrelated bacteria

Ready the workspace

Load all necessary libraries, packages, and data files

library(survival)
library(RColorBrewer)
library(ggplot2)
library(rlang)
library(agricolae)
library(coxme)
library(survminer)
library(olsrr)
library(tidyr)
library(dplyr)
library(Hmisc)
library(qqplotr)
library(ggthemes)
library(fabricatr)
library(gridExtra)
library(grid)
library(kableExtra)
library(sjPlot)
library(cowplot)
library(ggpubr)
library(patchwork)
library(SciViews)

rm(list=ls()) #removes all old variables

#setwd("/Volumes/GoogleDrive/My Drive/Research/Moria Chambers Fruit Flies/061322/Manuscript/")

Sm_Pr = read.csv("Sm_Pr_Figure1.csv")
Ef_Pr = read.csv("Ef_Pr_Figure2.csv")
Ef_bacterial_load = read.csv("Ef bacterial load_Figure3.csv")
Sm_bacterial_load = read.csv("Sm bacterial load_Figure3.csv")

Sm_Ps = read.csv("Sm_Ps_Figure4.csv")
Pr_varied_primary = read.csv("Varied primary Sm_Pr_Figure5.csv")
Ps_varied_primary = read.csv("Varied primary Sm_Ps_Figure5.csv")

Q1 Does the protection conferred by primary infection with S. marcescens or E. faecalis depend on the infectious dose of the secondary infection with P. rettgeri?

Study design:

  • Previous work: It has been previously published that chronic infection with S. marcescens and E. faecalis provide significant protection against secondary infection with P. rettgeri initiated by an infection dose of 23nL of 0.1 Abs600nm (~8000 CFU/fly) (Chambers et al. 2019)
  • We used infectious doses at least 10 fold above and below this dose, to assess how sensitive this protection is to dose.
  • Flies were first infected with S. marcescens or E. faecalis and then a week later were given P. rettgeri infections with a range of doses (0.001-5.0 Abs600nm) and survival monitored for at least three days.
  • Initially analyzed using Cox Proportional Hazards mixed effects model with date of injection, primary injection and secondary dose as fixed effects and random effect terms accounting for experimental structure (e.g. date of injection, housing of groups of flies in vials). However, this model violated the assumptions of Cox Proportional Hazards.
  • Therefore we switched to using log-rank pairwise comparisons to determine whether chronic infection significantly impacted survival of secondary infection at an individual dose.
    • LogRank comparisons were done on data from each date individually
    • If experiments on all dates showed higher survival in conditions with chronic infection, p-values were combined using Fisher’s test and only the final p-value is reported in the text.
    • If there was a variable effect of chronic infection, the number of individual experiments with significant effects is reported.
  • Conclusions
    • S. marcescens chronic infection provided significant protection during secondary infection with P. rettgeri at all infectious doses (p<0.001)
    • E. faecalis chronic infection with offered significant protection against three out of the four doses (P < 0.05). However, the one dose that did not show significant protection (Abs600nm = 0.1) had a borderline p-value = 0.06

Figure 1 - S. marcescens Chronic Infection

  • Mortality over 7 days following secondary infection with P. rettgeri in flies chronically infected with S. marcescens versus control flies not given a chronic infection.

Total Graph

#Sm_Pr$primary <- ordered(Sm_Pr$primary, levels = c("PBS", "S.m"))
table(Sm_Pr$primary, Sm_Pr$secondary_dose)
##      
##         0 0.01 0.1   1   2   3   5
##   PBS 115  118 206 283 263  80 165
##   S.m 114  119 196 275 238  81 159
#Sm_bacterial_load$primary <- ordered(Sm_bacterial_load$primary, levels = c("PBS", "S.m"))
#table(Sm_bacterial_load$primary, Sm_bacterial_load$secondary_dose)

#THE GRAPH
#png(file="Sm_Pr.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(Sm_Pr$day, Sm_Pr$status) ~ Sm_Pr$primary + Sm_Pr$secondary_dose), col=c("gray", "#D7EEF7", "#B5E8FC",  "#09B6FC", "#055df5", "#003694", "black"), lty=c(3,3,3,3, 3, 3, 3, 1,1,1,1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-PBS","PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0","PBS-P.rett 3.0","PBS-P.rett 5.0","S.m-PBS", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0", "S.m-P.rett 3.0", "S.m-P.rett 5.0"),col=c("gray", "#D7EEF7", "#B5E8FC",  "#09B6FC", "#055df5", "#003694", "black"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)

Appendix Figure 1A - Controls

Sm_Pr_PBS<-subset(Sm_Pr,secondary_dose=="0")
#table(Sm_Pr_PBS$primary)

#THE GRAPH
#png(file="AFigure1A.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Pr_PBS), col=c("gray"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("bottomleft",inset=0, bty="n" , c("double injection control (n=115)", "chronic infection only (n=114)"),col=c("gray"),lty=c(3, 1),lwd=5, cex=1)

#dev.off()

Figure 1 Panels

Figure 1A - 0.01

Sm_Pr_0.01<-subset(Sm_Pr,secondary_dose=="0.01")
#table(Sm_Pr_0.01$primary)

#THE GRAPH
#png(file="Figure1A.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Pr_0.01), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=118)", "chronic & secondary infection (n=119)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.7)

#dev.off()

Figure 1B - 0.1

Sm_Pr_0.1<-subset(Sm_Pr,secondary_dose=="0.1")
#table(Sm_Pr_0.1$primary)
#THE GRAPH
#png(file="Figure1B.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Pr_0.1), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=206)", "chronic & secondary infection (n=196)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.7)

#dev.off()

Figure 1c - 1.0

Sm_Pr_1.0<-subset(Sm_Pr,secondary_dose=="1")
#table(Sm_Pr_1.0$primary)
#THE GRAPH
#png(file="Figure1C.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Pr_1.0), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=283)", "chronic & secondary infection (n=275)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.7)

#dev.off()

Figure 1D - 2.0

Sm_Pr_2.0<-subset(Sm_Pr,secondary_dose=="2")
#table(Sm_Pr_2.0$primary)
#THE GRAPH
#png(file="Figure1D.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Pr_2.0), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=263)", "chronic & secondary infection (n=238)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.7)

#dev.off()

Figure 1E - 3.0

Sm_Pr_3.0<-subset(Sm_Pr,secondary_dose=="3")
#table(Sm_Pr_3.0$primary)
#THE GRAPH
#png(file="Figure1E.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Pr_3.0), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=80)", "chronic & secondary infection (n=81)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.7)

#dev.off()

Figure 1F - 5.0

Sm_Pr_5.0<-subset(Sm_Pr,secondary_dose=="5")
#table(Sm_Pr_5.0$primary)
#THE GRAPH
#png(file="Figure1F.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Pr_5.0), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=165)", "chronic & secondary infection (n=159)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.7)

#dev.off()

Figure 1

a <- ggdraw() +
  draw_image("Figure1A.png")

b <- ggdraw() +
  draw_image("Figure1B.png")

c <- ggdraw() +
  draw_image("Figure1C.png")

d <- ggdraw() +
  draw_image("Figure1D.png")

e <- ggdraw() +
  draw_image("Figure1E.png")

f <- ggdraw() +
  draw_image("Figure1F.png")

multiplot <- a + b + c + d + e + f
z <- multiplot + 
  plot_annotation(tag_levels = 'A')
ggsave("Fig.1.ChronicSm_Pr.pdf", z, device = "pdf")
z

Individual Blocks

Block 1

Sm_Pr_block1 <- subset(Sm_Pr, date == "8/2/2018")
table(Sm_Pr_block1$primary, Sm_Pr_block1$secondary_dose)
##      
##       0.01 0.1  1
##   PBS   17  18 18
##   S.m   15  15 15
Sm_Pr_0.01_1<-subset(Sm_Pr_block1, secondary_dose=="0.01")
#table(Sm_Pr_0.01_1$primary)
Sm_Pr_0.1_1<-subset(Sm_Pr_block1, secondary_dose=="0.1")
#table(Sm_Pr_0.1_1$primary)
Sm_Pr_1.0_1<-subset(Sm_Pr_block1, secondary_dose=="1")
#table(Sm_Pr_1.0_1$primary)
#Sm_Pr_2.0_1<-subset(Sm_Pr_block1, secondary_dose=="2")
#table(Sm_Pr_2.0_1$primary)
#Sm_Pr_3.0_1<-subset(Sm_Pr_block1, secondary_dose=="3")
#table(Sm_Pr_2.0_1$primary)
#Sm_Pr_5.0_1<-subset(Sm_Pr_block1, secondary_dose=="5")
#table(Sm_Pr_5.0_1$primary)

surv_Sm_Pr_0.01_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.01_1)
surv_Sm_Pr_0.01_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.01_1)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 17        8      5.3      1.37      3.19
## primary=S.m 15        3      5.7      1.28      3.19
## 
##  Chisq= 3.2  on 1 degrees of freedom, p= 0.07
surv_Sm_Pr_0.1_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_1)
surv_Sm_Pr_0.1_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_1)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 18       18     15.5     0.395      2.33
## primary=S.m 15       12     14.5     0.423      2.33
## 
##  Chisq= 2.3  on 1 degrees of freedom, p= 0.1
surv_Sm_Pr_1.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_1)
surv_Sm_Pr_1.0_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_1)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 18       18     15.8     0.301       5.3
## primary=S.m 15       13     15.2     0.314       5.3
## 
##  Chisq= 5.3  on 1 degrees of freedom, p= 0.02
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block1$day, Sm_Pr_block1$status) ~ Sm_Pr_block1$primary + Sm_Pr_block1$secondary_dose), col=c("#B5E8FC",  "#09B6FC", "#055df5"), lty=c(3, 3, 3, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0"),col=c("#B5E8FC", "#09B6FC", "#055df5"),lty=c(3, 3, 3, 1, 1, 1),lwd=5, cex=1)

Block 2

Sm_Pr_block2 <- subset(Sm_Pr, date == "8/6/2018")
table(Sm_Pr_block2$primary, Sm_Pr_block2$secondary_dose)
##      
##       0.01 0.1  1
##   PBS   19  27 12
##   S.m   24  23 15
Sm_Pr_0.01_2<-subset(Sm_Pr_block2, secondary_dose=="0.01")
#table(Sm_Pr_0.01_2$primary)
Sm_Pr_0.1_2<-subset(Sm_Pr_block2, secondary_dose=="0.1")
#table(Sm_Pr_0.1_2$primary)
Sm_Pr_1.0_2<-subset(Sm_Pr_block2, secondary_dose=="1")
#table(Sm_Pr_1.0_2$primary)
#Sm_Pr_2.0_1<-subset(Sm_Pr_block1, secondary_dose=="2")
#table(Sm_Pr_2.0_1$primary)

surv_Sm_Pr_0.01_2<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.01_2)
surv_Sm_Pr_0.01_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.01_2)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 19       13     8.24      2.74      5.57
## primary=S.m 24        9    13.76      1.64      5.57
## 
##  Chisq= 5.6  on 1 degrees of freedom, p= 0.02
surv_Sm_Pr_0.1_2<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_2)
surv_Sm_Pr_0.1_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_2)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27       20       13      3.72      10.5
## primary=S.m 23        8       15      3.24      10.5
## 
##  Chisq= 10.5  on 1 degrees of freedom, p= 0.001
surv_Sm_Pr_1.0_2<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_2)
surv_Sm_Pr_1.0_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_2)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 12       12     8.89     1.089      4.91
## primary=S.m 15       10    13.11     0.738      4.91
## 
##  Chisq= 4.9  on 1 degrees of freedom, p= 0.03
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block2$day, Sm_Pr_block2$status) ~ Sm_Pr_block2$primary + Sm_Pr_block2$secondary_dose), col=c("#B5E8FC",  "#09B6FC", "#055df5"), lty=c(3, 3, 3, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0"),col=c("#B5E8FC", "#09B6FC", "#055df5"),lty=c(3, 3, 3, 1, 1, 1),lwd=5, cex=1)

Block 3

Sm_Pr_block3 <- subset(Sm_Pr, date == "6/14/2019")
table(Sm_Pr_block3$primary, Sm_Pr_block3$secondary_dose)
##      
##       0.01 0.1  1  2
##   PBS   15  14 17 20
##   S.m   15  19 17 18
Sm_Pr_0.01_3<-subset(Sm_Pr_block3, secondary_dose=="0.01")
#table(Sm_Pr_0.01_3$primary)
Sm_Pr_0.1_3<-subset(Sm_Pr_block3, secondary_dose=="0.1")
#table(Sm_Pr_0.1_3$primary)
Sm_Pr_1.0_3<-subset(Sm_Pr_block3, secondary_dose=="1")
#table(Sm_Pr_1.0_3$primary)
Sm_Pr_2.0_3<-subset(Sm_Pr_block3, secondary_dose=="2")
#table(Sm_Pr_2.0_3$primary)

surv_Sm_Pr_0.01_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.01_3)
surv_Sm_Pr_0.01_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.01_3)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 15        9     6.24      1.22      3.41
## primary=S.m 15        4     6.76      1.13      3.41
## 
##  Chisq= 3.4  on 1 degrees of freedom, p= 0.06
surv_Sm_Pr_0.1_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_3)
surv_Sm_Pr_0.1_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_3)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 14        7     2.74      6.63      12.1
## primary=S.m 19        0     4.26      4.26      12.1
## 
##  Chisq= 12.1  on 1 degrees of freedom, p= 5e-04
surv_Sm_Pr_1.0_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_3)
surv_Sm_Pr_1.0_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_3)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 17       17      8.5      8.50        33
## primary=S.m 17        3     11.5      6.28        33
## 
##  Chisq= 33  on 1 degrees of freedom, p= 9e-09
surv_Sm_Pr_2.0_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_3)
surv_Sm_Pr_2.0_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_3)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 20       20     12.1      5.15      26.8
## primary=S.m 18        5     12.9      4.83      26.8
## 
##  Chisq= 26.8  on 1 degrees of freedom, p= 2e-07
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block3$day, Sm_Pr_block3$status) ~ Sm_Pr_block3$primary + Sm_Pr_block3$secondary_dose), col=c("#B5E8FC",  "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)

Block 4

Sm_Pr_block4 <- subset(Sm_Pr, date == "6/17/2019")
table(Sm_Pr_block4$primary, Sm_Pr_block4$secondary_dose)
##      
##       0.01 0.1  1  2
##   PBS   15  19 17 17
##   S.m   18  13 16 11
Sm_Pr_0.01_4<-subset(Sm_Pr_block4, secondary_dose=="0.01")
#table(Sm_Pr_0.01_4$primary)
Sm_Pr_0.1_4<-subset(Sm_Pr_block4, secondary_dose=="0.1")
#table(Sm_Pr_0.1_4$primary)
Sm_Pr_1.0_4<-subset(Sm_Pr_block4, secondary_dose=="1")
#table(Sm_Pr_1.0_4$primary)
Sm_Pr_2.0_4<-subset(Sm_Pr_block4, secondary_dose=="2")
#table(Sm_Pr_2.0_4$primary)

surv_Sm_Pr_0.01_4<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.01_4)
surv_Sm_Pr_0.01_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.01_4)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 15       10     7.59     0.767      1.76
## primary=S.m 18        8    10.41     0.559      1.76
## 
##  Chisq= 1.8  on 1 degrees of freedom, p= 0.2
surv_Sm_Pr_0.1_4<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_4)
surv_Sm_Pr_0.1_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_4)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 19        9     7.85     0.170     0.438
## primary=S.m 13        5     6.15     0.216     0.438
## 
##  Chisq= 0.4  on 1 degrees of freedom, p= 0.5
surv_Sm_Pr_1.0_4<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_4)
surv_Sm_Pr_1.0_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_4)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 17       16     9.76      4.00      15.7
## primary=S.m 16        6    12.24      3.18      15.7
## 
##  Chisq= 15.7  on 1 degrees of freedom, p= 7e-05
surv_Sm_Pr_2.0_4<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_4)
surv_Sm_Pr_2.0_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_4)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 17       17    12.14      1.94      16.7
## primary=S.m 11        5     9.86      2.39      16.7
## 
##  Chisq= 16.7  on 1 degrees of freedom, p= 4e-05
#surv_Sm_Pr_3.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_1)
#surv_Sm_Pr_3.0_1

#surv_Sm_Pr_5.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_1)
#surv_Sm_Pr_5.0_1

#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block4$day, Sm_Pr_block4$status) ~ Sm_Pr_block4$primary + Sm_Pr_block4$secondary_dose), col=c("#B5E8FC",  "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)

Block 5

Sm_Pr_block5 <- subset(Sm_Pr, date == "6/19/2019")
table(Sm_Pr_block5$primary, Sm_Pr_block5$secondary_dose)
##      
##       0.01 0.1  1  2
##   PBS   15   5  9 19
##   S.m   11  17 14 15
Sm_Pr_0.01_5<-subset(Sm_Pr_block5, secondary_dose=="0.01")
#table(Sm_Pr_0.01_5$primary)
Sm_Pr_0.1_5<-subset(Sm_Pr_block5, secondary_dose=="0.1")
#table(Sm_Pr_0.1_5$primary)
Sm_Pr_1.0_5<-subset(Sm_Pr_block5, secondary_dose=="1")
#table(Sm_Pr_1.0_5$primary)
Sm_Pr_2.0_5<-subset(Sm_Pr_block5, secondary_dose=="2")
#table(Sm_Pr_2.0_5$primary)

surv_Sm_Pr_0.01_5<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.01_5)
surv_Sm_Pr_0.01_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.01_5)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 15       10     7.83     0.604      1.61
## primary=S.m 11        5     7.17     0.659      1.61
## 
##  Chisq= 1.6  on 1 degrees of freedom, p= 0.2
surv_Sm_Pr_0.1_5<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_5)
surv_Sm_Pr_0.1_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_5)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS  5        5     1.21     11.90      17.1
## primary=S.m 17        5     8.79      1.64      17.1
## 
##  Chisq= 17.1  on 1 degrees of freedom, p= 4e-05
surv_Sm_Pr_1.0_5<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_5)
surv_Sm_Pr_1.0_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_5)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS  9        9     5.29      2.61       8.8
## primary=S.m 14        7    10.71      1.29       8.8
## 
##  Chisq= 8.8  on 1 degrees of freedom, p= 0.003
surv_Sm_Pr_2.0_5<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_5)
surv_Sm_Pr_2.0_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_5)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 19       18     11.8      3.27        15
## primary=S.m 15        9     15.2      2.54        15
## 
##  Chisq= 15  on 1 degrees of freedom, p= 1e-04
#surv_Sm_Pr_3.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_1)
#surv_Sm_Pr_3.0_1

#surv_Sm_Pr_5.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_1)
#surv_Sm_Pr_5.0_1

#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block5$day, Sm_Pr_block5$status) ~ Sm_Pr_block5$primary + Sm_Pr_block5$secondary_dose), col=c("#B5E8FC",  "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)

Block 6

Sm_Pr_block6 <- subset(Sm_Pr, date == "6/25/2019")
table(Sm_Pr_block6$primary, Sm_Pr_block6$secondary_dose)
##      
##       0.01 0.1  1  2
##   PBS   17  19 19 20
##   S.m   20  15 17 16
Sm_Pr_0.01_6<-subset(Sm_Pr_block6, secondary_dose=="0.01")
#table(Sm_Pr_0.01_6$primary)
Sm_Pr_0.1_6<-subset(Sm_Pr_block6, secondary_dose=="0.1")
#table(Sm_Pr_0.1_6$primary)
Sm_Pr_1.0_6<-subset(Sm_Pr_block6, secondary_dose=="1")
#table(Sm_Pr_1.0_6$primary)
Sm_Pr_2.0_6<-subset(Sm_Pr_block6, secondary_dose=="2")
#table(Sm_Pr_2.0_6$primary)

surv_Sm_Pr_0.01_6<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.01_6)
surv_Sm_Pr_0.01_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.01_6)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 17       12     7.11      3.36      7.48
## primary=S.m 20        5     9.89      2.42      7.48
## 
##  Chisq= 7.5  on 1 degrees of freedom, p= 0.006
surv_Sm_Pr_0.1_6<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_6)
surv_Sm_Pr_0.1_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_6)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 19       15     8.52      4.92      13.2
## primary=S.m 15        2     8.48      4.95      13.2
## 
##  Chisq= 13.2  on 1 degrees of freedom, p= 3e-04
surv_Sm_Pr_1.0_6<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_6)
surv_Sm_Pr_1.0_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_6)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 19       17     9.99      4.92        18
## primary=S.m 17        5    12.01      4.09        18
## 
##  Chisq= 18  on 1 degrees of freedom, p= 2e-05
surv_Sm_Pr_2.0_6<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_6)
surv_Sm_Pr_2.0_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_6)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 20       20     11.9      5.57      27.8
## primary=S.m 16        6     14.1      4.68      27.8
## 
##  Chisq= 27.8  on 1 degrees of freedom, p= 1e-07
#surv_Sm_Pr_3.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_1)
#surv_Sm_Pr_3.0_1

#surv_Sm_Pr_5.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_1)
#surv_Sm_Pr_5.0_1

#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block6$day, Sm_Pr_block6$status) ~ Sm_Pr_block6$primary + Sm_Pr_block6$secondary_dose), col=c("#B5E8FC",  "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)

Block 7

Sm_Pr_block7 <- subset(Sm_Pr, date == "6/27/2019")
table(Sm_Pr_block7$primary, Sm_Pr_block7$secondary_dose)
##      
##       0.01 0.1  1  2
##   PBS   20  17 20 19
##   S.m   16  14 14 14
Sm_Pr_0.01_7<-subset(Sm_Pr_block7, secondary_dose=="0.01")
#table(Sm_Pr_0.01_7$primary)
Sm_Pr_0.1_7<-subset(Sm_Pr_block7, secondary_dose=="0.1")
#table(Sm_Pr_0.1_7$primary)
Sm_Pr_1.0_7<-subset(Sm_Pr_block7, secondary_dose=="1")
#table(Sm_Pr_1.0_7$primary)
Sm_Pr_2.0_7<-subset(Sm_Pr_block7, secondary_dose=="2")
#table(Sm_Pr_2.0_7$primary)

surv_Sm_Pr_0.01_7<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.01_7)
surv_Sm_Pr_0.01_7
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.01_7)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 20       11    10.11    0.0779     0.206
## primary=S.m 16        9     9.89    0.0797     0.206
## 
##  Chisq= 0.2  on 1 degrees of freedom, p= 0.7
surv_Sm_Pr_0.1_7<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_7)
surv_Sm_Pr_0.1_7
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_7)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 17       14     9.38      2.27      6.94
## primary=S.m 14        6    10.62      2.01      6.94
## 
##  Chisq= 6.9  on 1 degrees of freedom, p= 0.008
surv_Sm_Pr_1.0_7<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_7)
surv_Sm_Pr_1.0_7
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_7)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 20       20     13.5      3.09      22.5
## primary=S.m 14        7     13.5      3.11      22.5
## 
##  Chisq= 22.5  on 1 degrees of freedom, p= 2e-06
surv_Sm_Pr_2.0_7<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_7)
surv_Sm_Pr_2.0_7
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_7)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 19       18     11.2      4.18      17.4
## primary=S.m 14        9     15.8      2.95      17.4
## 
##  Chisq= 17.4  on 1 degrees of freedom, p= 3e-05
#surv_Sm_Pr_3.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_1)
#surv_Sm_Pr_3.0_1

#surv_Sm_Pr_5.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_1)
#surv_Sm_Pr_5.0_1

#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block7$day, Sm_Pr_block7$status) ~ Sm_Pr_block7$primary + Sm_Pr_block7$secondary_dose), col=c("#B5E8FC",  "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)

Block 8

Sm_Pr_block8 <- subset(Sm_Pr, date == "1/30/2020")
table(Sm_Pr_block8$primary, Sm_Pr_block8$secondary_dose)
##      
##        0 0.1  1  2  5
##   PBS 18  29 27 26 28
##   S.m 20  23 24 25 25
Sm_Pr_0.1_8<-subset(Sm_Pr_block8, secondary_dose=="0.1")
#table(Sm_Pr_0.1_8$primary)
Sm_Pr_1.0_8<-subset(Sm_Pr_block8, secondary_dose=="1")
#table(Sm_Pr_1.0_8$primary)
Sm_Pr_2.0_8<-subset(Sm_Pr_block8, secondary_dose=="2")
#table(Sm_Pr_2.0_8$primary)
#Sm_Pr_3.0_8<-subset(Sm_Pr_block8, secondary_dose=="3")
#table(Sm_Pr_3.0_8$primary)
Sm_Pr_5.0_8<-subset(Sm_Pr_block8, secondary_dose=="5")
#table(Sm_Pr_5.0_8$primary)


surv_Sm_Pr_0.1_8<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_8)
surv_Sm_Pr_0.1_8
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_8)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 29       29     15.7     11.33      26.8
## primary=S.m 23       15     28.3      6.27      26.8
## 
##  Chisq= 26.8  on 1 degrees of freedom, p= 2e-07
surv_Sm_Pr_1.0_8<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_8)
surv_Sm_Pr_1.0_8
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_8)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27       27     16.9      5.97      33.4
## primary=S.m 24       19     29.1      3.48      33.4
## 
##  Chisq= 33.4  on 1 degrees of freedom, p= 8e-09
surv_Sm_Pr_2.0_8<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_8)
surv_Sm_Pr_2.0_8
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_8)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 26       26     15.3      7.49      36.4
## primary=S.m 25       20     30.7      3.73      36.4
## 
##  Chisq= 36.4  on 1 degrees of freedom, p= 2e-09
#surv_Sm_Pr_3.0_8<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_8)
#surv_Sm_Pr_3.0_8

surv_Sm_Pr_5.0_8<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_8)
surv_Sm_Pr_5.0_8
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_5.0_8)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28       28     15.3     10.49      48.2
## primary=S.m 25       16     28.7      5.61      48.2
## 
##  Chisq= 48.2  on 1 degrees of freedom, p= 4e-12
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block8$day, Sm_Pr_block8$status) ~ Sm_Pr_block8$primary + Sm_Pr_block8$secondary_dose), col=c("gray", "#B5E8FC",  "#09B6FC", "#055df5", "black"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-PBS", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "PBS-P.rett 5.0", "S.m-PBS", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0", "S.m-P.rett 5.0"),col=c("gray", "#B5E8FC",  "#09B6FC", "#055df5", "black"),lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=5, cex=1)

Block 9

Sm_Pr_block9 <- subset(Sm_Pr, date == "2/13/2020")
table(Sm_Pr_block9$primary, Sm_Pr_block9$secondary_dose)
##      
##        0 0.1  1  2  5
##   PBS 18  28 31 30 30
##   S.m 18  28 28 28 28
Sm_Pr_0.1_9<-subset(Sm_Pr_block9, secondary_dose=="0.1")
#table(Sm_Pr_0.1_9$primary)
Sm_Pr_1.0_9<-subset(Sm_Pr_block9, secondary_dose=="1")
#table(Sm_Pr_1.0_9$primary)
Sm_Pr_2.0_9<-subset(Sm_Pr_block9, secondary_dose=="2")
#table(Sm_Pr_2.0_9$primary)
#Sm_Pr_3.0_9<-subset(Sm_Pr_block9, secondary_dose=="3")
#table(Sm_Pr_3.0_9$primary)
Sm_Pr_5.0_9<-subset(Sm_Pr_block9, secondary_dose=="5")
#table(Sm_Pr_5.0_9$primary)


surv_Sm_Pr_0.1_9<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_9)
surv_Sm_Pr_0.1_9
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_9)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28       26     13.9     10.64      21.8
## primary=S.m 28       10     22.1      6.66      21.8
## 
##  Chisq= 21.8  on 1 degrees of freedom, p= 3e-06
surv_Sm_Pr_1.0_9<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_9)
surv_Sm_Pr_1.0_9
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_9)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 31       31     18.9      7.72        41
## primary=S.m 28       15     27.1      5.39        41
## 
##  Chisq= 41  on 1 degrees of freedom, p= 2e-10
surv_Sm_Pr_2.0_9<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_9)
surv_Sm_Pr_2.0_9
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_9)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30       30     17.1      9.80      46.3
## primary=S.m 28        8     20.9      7.99      46.3
## 
##  Chisq= 46.3  on 1 degrees of freedom, p= 1e-11
#surv_Sm_Pr_3.0_9<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_9)
#surv_Sm_Pr_3.0_9

surv_Sm_Pr_5.0_9<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_9)
surv_Sm_Pr_5.0_9
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_5.0_9)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30       30     18.2      7.69      35.3
## primary=S.m 28       21     32.8      4.26      35.3
## 
##  Chisq= 35.3  on 1 degrees of freedom, p= 3e-09
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block9$day, Sm_Pr_block9$status) ~ Sm_Pr_block9$primary + Sm_Pr_block9$secondary_dose), col=c("gray", "#B5E8FC",  "#09B6FC", "#055df5", "black"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-PBS", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "PBS-P.rett 5.0", "S.m-PBS", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0", "S.m-P.rett 5.0"),col=c("gray", "#B5E8FC",  "#09B6FC", "#055df5", "black"),lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=5, cex=1)

Block 10

Sm_Pr_block10 <- subset(Sm_Pr, date == "9/27/2020")
table(Sm_Pr_block10$primary, Sm_Pr_block10$secondary_dose)
##      
##        0  1  2  3  5
##   PBS 10 27 26 20 20
##   S.m 10 28 27 28 27
#Sm_Pr_0.1_10<-subset(Sm_Pr_block10, secondary_dose=="0.1")
#table(Sm_Pr_0.1_10$primary)
Sm_Pr_1.0_10<-subset(Sm_Pr_block10, secondary_dose=="1")
#table(Sm_Pr_1.0_10$primary)
Sm_Pr_2.0_10<-subset(Sm_Pr_block10, secondary_dose=="2")
#table(Sm_Pr_2.0_10$primary)
Sm_Pr_3.0_10<-subset(Sm_Pr_block10, secondary_dose=="3")
#table(Sm_Pr_3.0_10$primary)
Sm_Pr_5.0_10<-subset(Sm_Pr_block10, secondary_dose=="5")
#table(Sm_Pr_5.0_10$primary)


#surv_Sm_Pr_0.1_10<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_10)
#surv_Sm_Pr_0.1_10


surv_Sm_Pr_1.0_10<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_10)
surv_Sm_Pr_1.0_10
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_10)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27       27     12.7     16.16      53.6
## primary=S.m 28        7     21.3      9.62      53.6
## 
##  Chisq= 53.6  on 1 degrees of freedom, p= 2e-13
surv_Sm_Pr_2.0_10<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_10)
surv_Sm_Pr_2.0_10
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_10)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 26       26     14.7      8.65      38.4
## primary=S.m 27        9     20.3      6.28      38.4
## 
##  Chisq= 38.4  on 1 degrees of freedom, p= 6e-10
surv_Sm_Pr_3.0_10<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_10)
surv_Sm_Pr_3.0_10
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_3.0_10)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 20       20     10.4      8.82      30.9
## primary=S.m 28       10     19.6      4.69      30.9
## 
##  Chisq= 30.9  on 1 degrees of freedom, p= 3e-08
surv_Sm_Pr_5.0_10<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_10)
surv_Sm_Pr_5.0_10
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_5.0_10)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 20       20     12.8      4.10      19.3
## primary=S.m 27       16     23.2      2.25      19.3
## 
##  Chisq= 19.3  on 1 degrees of freedom, p= 1e-05
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block10$day, Sm_Pr_block10$status) ~ Sm_Pr_block10$primary + Sm_Pr_block10$secondary_dose), col=c("gray", "#09B6FC", "#055df5", "#003694", "black"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-PBS", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "PBS-P.rett 3.0", "PBS-P.rett 5.0", "S.m-PBS", "S.m-P.rett 1.0", "S.m-P.rett 2.0", "S.m-P.rett 3.0", "S.m-P.rett 5.0"),col=c("gray", "#09B6FC", "#055df5", "#003694", "black"),lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=5, cex=1)

Block 11

Sm_Pr_block11 <- subset(Sm_Pr, date == "10/6/2020")
table(Sm_Pr_block11$primary, Sm_Pr_block11$secondary_dose)
##      
##        0  1  2  3  5
##   PBS 20 28 28 30 29
##   S.m 16 28 27 26 27
#Sm_Pr_0.1_11<-subset(Sm_Pr_block11, secondary_dose=="0.1")
#table(Sm_Pr_0.1_11$primary)
Sm_Pr_1.0_11<-subset(Sm_Pr_block11, secondary_dose=="1")
#table(Sm_Pr_1.0_11$primary)
Sm_Pr_2.0_11<-subset(Sm_Pr_block11, secondary_dose=="2")
#table(Sm_Pr_2.0_11$primary)
Sm_Pr_3.0_11<-subset(Sm_Pr_block11, secondary_dose=="3")
#table(Sm_Pr_3.0_11$primary)
Sm_Pr_5.0_11<-subset(Sm_Pr_block11, secondary_dose=="5")
#table(Sm_Pr_5.0_11$primary)


#surv_Sm_Pr_0.1_11<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_11)
#surv_Sm_Pr_0.1_11


surv_Sm_Pr_1.0_11<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_11)
surv_Sm_Pr_1.0_11
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_11)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28       24     14.2      6.85        19
## primary=S.m 28       12     21.8      4.44        19
## 
##  Chisq= 19  on 1 degrees of freedom, p= 1e-05
surv_Sm_Pr_2.0_11<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_11)
surv_Sm_Pr_2.0_11
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_11)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28       28     14.4     12.93      51.5
## primary=S.m 27        9     22.6      8.21      51.5
## 
##  Chisq= 51.5  on 1 degrees of freedom, p= 7e-13
surv_Sm_Pr_3.0_11<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_11)
surv_Sm_Pr_3.0_11
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_3.0_11)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30       30     16.1     12.07        55
## primary=S.m 26       14     27.9      6.95        55
## 
##  Chisq= 55  on 1 degrees of freedom, p= 1e-13
surv_Sm_Pr_5.0_11<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_11)
surv_Sm_Pr_5.0_11
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_5.0_11)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 29       29     16.1      10.4      47.6
## primary=S.m 27       11     23.9       7.0      47.6
## 
##  Chisq= 47.6  on 1 degrees of freedom, p= 5e-12
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block11$day, Sm_Pr_block11$status) ~ Sm_Pr_block11$primary + Sm_Pr_block11$secondary_dose), col=c("gray", "#09B6FC", "#055df5", "#003694", "black"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-PBS", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "PBS-P.rett 3.0", "PBS-P.rett 5.0", "S.m-PBS", "S.m-P.rett 1.0", "S.m-P.rett 2.0", "S.m-P.rett 3.0", "S.m-P.rett 5.0"),col=c("gray", "#09B6FC", "#055df5", "#003694", "black"),lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=5, cex=1)

Block 12

Sm_Pr_block12 <- subset(Sm_Pr, date == "11/1/2020")
table(Sm_Pr_block12$primary, Sm_Pr_block12$secondary_dose)
##      
##        0  1  2  3  5
##   PBS 20 30 30 30 30
##   S.m 20 30 30 27 25
#Sm_Pr_0.1_12<-subset(Sm_Pr_block12, secondary_dose=="0.1")
#table(Sm_Pr_0.1_12$primary)
Sm_Pr_1.0_12<-subset(Sm_Pr_block12, secondary_dose=="1")
#table(Sm_Pr_1.0_12$primary)
Sm_Pr_2.0_12<-subset(Sm_Pr_block12, secondary_dose=="2")
#table(Sm_Pr_2.0_12$primary)
Sm_Pr_3.0_12<-subset(Sm_Pr_block12, secondary_dose=="3")
#table(Sm_Pr_3.0_12$primary)
Sm_Pr_5.0_12<-subset(Sm_Pr_block12, secondary_dose=="5")
#table(Sm_Pr_5.0_12$primary)


#surv_Sm_Pr_0.1_12<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_12)
#surv_Sm_Pr_0.1_12


surv_Sm_Pr_1.0_12<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_12)
surv_Sm_Pr_1.0_12
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_12)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30       28     13.6     15.17      33.2
## primary=S.m 30       13     27.4      7.55      33.2
## 
##  Chisq= 33.2  on 1 degrees of freedom, p= 8e-09
surv_Sm_Pr_2.0_12<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_12)
surv_Sm_Pr_2.0_12
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_12)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30       30     13.9      18.5      58.6
## primary=S.m 30       10     26.1       9.9      58.6
## 
##  Chisq= 58.6  on 1 degrees of freedom, p= 2e-14
surv_Sm_Pr_3.0_12<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_12)
surv_Sm_Pr_3.0_12
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_3.0_12)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30       30     15.6     13.39        50
## primary=S.m 27       21     35.4      5.88        50
## 
##  Chisq= 50  on 1 degrees of freedom, p= 2e-12
surv_Sm_Pr_5.0_12<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_12)
surv_Sm_Pr_5.0_12
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_5.0_12)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30       30     20.2      4.78      31.5
## primary=S.m 25       17     26.8      3.59      31.5
## 
##  Chisq= 31.5  on 1 degrees of freedom, p= 2e-08
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block12$day, Sm_Pr_block12$status) ~ Sm_Pr_block12$primary + Sm_Pr_block12$secondary_dose), col=c("gray", "#09B6FC", "#055df5", "#003694", "black"), lty=c(3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-PBS", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "PBS-P.rett 3.0", "PBS-P.rett 5.0", "S.m-PBS", "S.m-P.rett 1.0", "S.m-P.rett 2.0", "S.m-P.rett 3.0", "S.m-P.rett 5.0"),col=c("gray", "#09B6FC", "#055df5", "#003694", "black"),lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=5, cex=1)

Block 13

Sm_Pr_block13 <- subset(Sm_Pr, date == "6/16/2021")
table(Sm_Pr_block13$primary, Sm_Pr_block13$secondary_dose)
##      
##        0 0.1  1  2  5
##   PBS 29  30 28 28 28
##   S.m 30  29 29 27 27
Sm_Pr_0.1_13<-subset(Sm_Pr_block13, secondary_dose=="0.1")
#table(Sm_Pr_0.1_13$primary)
Sm_Pr_1.0_13<-subset(Sm_Pr_block13, secondary_dose=="1")
#table(Sm_Pr_1.0_13$primary)
Sm_Pr_2.0_13<-subset(Sm_Pr_block13, secondary_dose=="2")
#table(Sm_Pr_2.0_13$primary)
#Sm_Pr_3.0_13<-subset(Sm_Pr_block13, secondary_dose=="3")
#table(Sm_Pr_2.0_13$primary)
Sm_Pr_5.0_13<-subset(Sm_Pr_block13, secondary_dose=="5")
#table(Sm_Pr_2.0_13$primary)


surv_Sm_Pr_0.1_13<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_13)
surv_Sm_Pr_0.1_13
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_13)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30       28     14.5     12.56      39.3
## primary=S.m 29        6     19.5      9.34      39.3
## 
##  Chisq= 39.3  on 1 degrees of freedom, p= 4e-10
surv_Sm_Pr_1.0_13<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_13)
surv_Sm_Pr_1.0_13
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_13)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28       28     14.3      13.2      51.4
## primary=S.m 29        4     17.7      10.6      51.4
## 
##  Chisq= 51.4  on 1 degrees of freedom, p= 7e-13
surv_Sm_Pr_2.0_13<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_13)
surv_Sm_Pr_2.0_13
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_13)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28       28     15.4     10.42      45.1
## primary=S.m 27        5     17.6      9.06      45.1
## 
##  Chisq= 45.1  on 1 degrees of freedom, p= 2e-11
#surv_Sm_Pr_3.0_13<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_13)
#surv_Sm_Pr_3.0_13

surv_Sm_Pr_5.0_13<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_13)
surv_Sm_Pr_5.0_13
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_5.0_13)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28       28     15.4     10.42      45.1
## primary=S.m 27        4     16.6      9.61      45.1
## 
##  Chisq= 45.1  on 1 degrees of freedom, p= 2e-11
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block13$day, Sm_Pr_block13$status) ~ Sm_Pr_block13$primary + Sm_Pr_block13$secondary_dose), col=c("gray", "#09B6FC", "#055df5", "#003694", "black"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-PBS", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "PBS-P.rett 3.0", "PBS-P.rett 5.0", "S.m-PBS", "S.m-P.rett 1.0", "S.m-P.rett 2.0", "S.m-P.rett 3.0", "S.m-P.rett 5.0"),col=c("gray", "#09B6FC", "#055df5", "#003694", "black"),lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=5, cex=1)

Statistical tests

# Grab p-values from each individual logrank test (results shown above in each block)
pvalue_0.01_1<-pchisq(surv_Sm_Pr_0.01_1$chisq, df=1, lower.tail=F)
pvalue_0.1_1<-pchisq(surv_Sm_Pr_0.1_1$chisq, df=1, lower.tail=F)
pvalue_1.0_1<-pchisq(surv_Sm_Pr_1.0_1$chisq, df=1, lower.tail=F)

pvalue_0.01_2<-pchisq(surv_Sm_Pr_0.01_2$chisq, df=1, lower.tail=F)
pvalue_0.1_2<-pchisq(surv_Sm_Pr_0.1_2$chisq, df=1, lower.tail=F)
pvalue_1.0_2<-pchisq(surv_Sm_Pr_1.0_2$chisq, df=1, lower.tail=F)

pvalue_0.01_3<-pchisq(surv_Sm_Pr_0.01_3$chisq, df=1, lower.tail=F)
pvalue_0.1_3<-pchisq(surv_Sm_Pr_0.1_3$chisq, df=1, lower.tail=F)
pvalue_1.0_3<-pchisq(surv_Sm_Pr_1.0_3$chisq, df=1, lower.tail=F)
pvalue_2.0_3<-pchisq(surv_Sm_Pr_2.0_3$chisq, df=1, lower.tail=F)

pvalue_0.01_4<-pchisq(surv_Sm_Pr_0.01_4$chisq, df=1, lower.tail=F)
pvalue_0.1_4<-pchisq(surv_Sm_Pr_0.1_4$chisq, df=1, lower.tail=F)
pvalue_1.0_4<-pchisq(surv_Sm_Pr_1.0_4$chisq, df=1, lower.tail=F)
pvalue_2.0_4<-pchisq(surv_Sm_Pr_2.0_4$chisq, df=1, lower.tail=F)

pvalue_0.01_5<-pchisq(surv_Sm_Pr_0.01_5$chisq, df=1, lower.tail=F)
pvalue_0.1_5<-pchisq(surv_Sm_Pr_0.1_5$chisq, df=1, lower.tail=F)
pvalue_1.0_5<-pchisq(surv_Sm_Pr_1.0_5$chisq, df=1, lower.tail=F)
pvalue_2.0_5<-pchisq(surv_Sm_Pr_2.0_5$chisq, df=1, lower.tail=F)

pvalue_0.01_6<-pchisq(surv_Sm_Pr_0.01_6$chisq, df=1, lower.tail=F)
pvalue_0.1_6<-pchisq(surv_Sm_Pr_0.1_6$chisq, df=1, lower.tail=F)
pvalue_1.0_6<-pchisq(surv_Sm_Pr_1.0_6$chisq, df=1, lower.tail=F)
pvalue_2.0_6<-pchisq(surv_Sm_Pr_2.0_6$chisq, df=1, lower.tail=F)

pvalue_0.01_7<-pchisq(surv_Sm_Pr_0.01_7$chisq, df=1, lower.tail=F)
pvalue_0.1_7<-pchisq(surv_Sm_Pr_0.1_7$chisq, df=1, lower.tail=F)
pvalue_1.0_7<-pchisq(surv_Sm_Pr_1.0_7$chisq, df=1, lower.tail=F)
pvalue_2.0_7<-pchisq(surv_Sm_Pr_2.0_7$chisq, df=1, lower.tail=F)


pvalue_0.1_8<-pchisq(surv_Sm_Pr_0.1_8$chisq, df=1, lower.tail=F)
pvalue_1.0_8<-pchisq(surv_Sm_Pr_1.0_8$chisq, df=1, lower.tail=F)
pvalue_2.0_8<-pchisq(surv_Sm_Pr_2.0_8$chisq, df=1, lower.tail=F)
pvalue_5.0_8<-pchisq(surv_Sm_Pr_5.0_8$chisq, df=1, lower.tail=F)


pvalue_0.1_9<-pchisq(surv_Sm_Pr_0.1_9$chisq, df=1, lower.tail=F)
pvalue_1.0_9<-pchisq(surv_Sm_Pr_1.0_9$chisq, df=1, lower.tail=F)
pvalue_2.0_9<-pchisq(surv_Sm_Pr_2.0_9$chisq, df=1, lower.tail=F)
pvalue_5.0_9<-pchisq(surv_Sm_Pr_5.0_9$chisq, df=1, lower.tail=F)

pvalue_1.0_10<-pchisq(surv_Sm_Pr_1.0_10$chisq, df=1, lower.tail=F)
pvalue_2.0_10<-pchisq(surv_Sm_Pr_2.0_10$chisq, df=1, lower.tail=F)
pvalue_3.0_10<-pchisq(surv_Sm_Pr_3.0_10$chisq, df=1, lower.tail=F)
pvalue_5.0_10<-pchisq(surv_Sm_Pr_5.0_10$chisq, df=1, lower.tail=F)

pvalue_1.0_11<-pchisq(surv_Sm_Pr_1.0_11$chisq, df=1, lower.tail=F)
pvalue_2.0_11<-pchisq(surv_Sm_Pr_2.0_11$chisq, df=1, lower.tail=F)
pvalue_3.0_11<-pchisq(surv_Sm_Pr_3.0_11$chisq, df=1, lower.tail=F)
pvalue_5.0_11<-pchisq(surv_Sm_Pr_5.0_11$chisq, df=1, lower.tail=F)

pvalue_1.0_12<-pchisq(surv_Sm_Pr_1.0_12$chisq, df=1, lower.tail=F)
pvalue_2.0_12<-pchisq(surv_Sm_Pr_2.0_12$chisq, df=1, lower.tail=F)
pvalue_3.0_12<-pchisq(surv_Sm_Pr_3.0_12$chisq, df=1, lower.tail=F)
pvalue_5.0_12<-pchisq(surv_Sm_Pr_5.0_12$chisq, df=1, lower.tail=F)

pvalue_0.1_13<-pchisq(surv_Sm_Pr_0.1_13$chisq, df=1, lower.tail=F)
pvalue_1.0_13<-pchisq(surv_Sm_Pr_1.0_13$chisq, df=1, lower.tail=F)
pvalue_2.0_13<-pchisq(surv_Sm_Pr_2.0_13$chisq, df=1, lower.tail=F)
pvalue_5.0_13<-pchisq(surv_Sm_Pr_5.0_13$chisq, df=1, lower.tail=F)
# Use all p-values from a given condition to get Fisher's combined test p-value
pvalues_0.01 <- c(pvalue_0.01_1,pvalue_0.01_2,pvalue_0.01_3,pvalue_0.01_4,pvalue_0.01_5,pvalue_0.01_6,pvalue_0.01_7)
test.stat_0.01 <- -2*sum(log(pvalues_0.01))
pfinal_0.01<-pchisq(test.stat_0.01, df=2*length(pvalues_0.01), lower.tail=F)

pvalues_0.1 <- c(pvalue_0.1_1,pvalue_0.1_2,pvalue_0.1_3,pvalue_0.1_4,pvalue_0.1_5,pvalue_0.1_6,pvalue_0.1_7, pvalue_0.1_8, pvalue_0.1_9, pvalue_0.1_13)
test.stat_0.1 <- -2*sum(log(pvalues_0.1))
pfinal_0.1<-pchisq(test.stat_0.1, df=2*length(pvalues_0.1), lower.tail=F)

pvalues_1.0 <- c(pvalue_1.0_1,pvalue_1.0_2,pvalue_1.0_3,pvalue_1.0_4,pvalue_1.0_5,pvalue_1.0_6,pvalue_1.0_7, pvalue_1.0_8, pvalue_1.0_9, pvalue_1.0_10, pvalue_1.0_11, pvalue_1.0_12, pvalue_1.0_13)
test.stat_1.0 <- -2*sum(log(pvalues_1.0))
pfinal_1.0<-pchisq(test.stat_1.0, df=2*length(pvalues_1.0), lower.tail=F)

pvalues_2.0 <- c(pvalue_2.0_3,pvalue_2.0_4,pvalue_2.0_5,pvalue_2.0_6,pvalue_2.0_7, pvalue_2.0_8, pvalue_2.0_9, pvalue_2.0_10, pvalue_2.0_11, pvalue_2.0_12, pvalue_2.0_13)
test.stat_2.0 <- -2*sum(log(pvalues_2.0))
pfinal_2.0<-pchisq(test.stat_2.0, df=2*length(pvalues_2.0), lower.tail=F)

pvalues_3.0 <- c(pvalue_3.0_10, pvalue_3.0_11, pvalue_3.0_12)
test.stat_3.0 <- -2*sum(log(pvalues_3.0))
pfinal_3.0<-pchisq(test.stat_3.0, df=2*length(pvalues_3.0), lower.tail=F)

pvalues_5.0 <- c(pvalue_5.0_8, pvalue_5.0_9, pvalue_5.0_10, pvalue_5.0_11, pvalue_5.0_12, pvalue_5.0_13)
test.stat_5.0 <- -2*sum(log(pvalues_5.0))
pfinal_5.0<-pchisq(test.stat_5.0, df=2*length(pvalues_5.0), lower.tail=F)
# Generate Table of final results
pvalues_0.01_table <- c(pvalue_0.01_1,pvalue_0.01_2,pvalue_0.01_3,pvalue_0.01_4,pvalue_0.01_5,pvalue_0.01_6,pvalue_0.01_7, NA, NA, NA, NA, NA, NA)
p0.01<-scales::pvalue(pvalues_0.01_table)
p0.01 <- p0.01 %>% replace_na("--")

pvalues_0.1_table <- c(pvalue_0.1_1,pvalue_0.1_2,pvalue_0.1_3,pvalue_0.1_4,pvalue_0.1_5,pvalue_0.1_6,pvalue_0.1_7, pvalue_0.1_8, pvalue_0.1_9,NA, NA, NA, pvalue_0.1_13)
p0.1<-scales::pvalue(pvalues_0.1_table)
p0.1<- p0.1 %>% replace_na("--")

pvalues_1.0_table <- c(pvalue_1.0_1,pvalue_1.0_2,pvalue_1.0_3,pvalue_1.0_4,pvalue_1.0_5,pvalue_1.0_6,pvalue_1.0_7, pvalue_1.0_8, pvalue_1.0_9, pvalue_1.0_10, pvalue_1.0_11, pvalue_1.0_12, pvalue_1.0_13)
p1.0<-scales::pvalue(pvalues_1.0_table)

pvalues_2.0_table <- c(NA, NA, pvalue_2.0_3, pvalue_2.0_4, pvalue_2.0_5, pvalue_2.0_6, pvalue_2.0_7, pvalue_2.0_8, pvalue_2.0_9, pvalue_2.0_10, pvalue_2.0_11, pvalue_2.0_12, pvalue_2.0_13)
p2.0<-scales::pvalue(pvalues_2.0_table)
p2.0<-p2.0 %>% replace_na("--")

pvalues_3.0_table <- c(NA, NA, NA, NA, NA, NA, NA, NA, NA, pvalue_3.0_10, pvalue_3.0_11, pvalue_3.0_12, NA)
p3.0<-scales::pvalue(pvalues_3.0_table)
p3.0<-p3.0 %>% replace_na("--")

pvalues_5.0_table <- c(NA, NA, NA, NA, NA, NA, NA, pvalue_5.0_8, pvalue_5.0_9, pvalue_5.0_10, pvalue_5.0_11, pvalue_5.0_12, pvalue_5.0_13)
p5.0<-scales::pvalue(pvalues_5.0_table)
p5.0<-p5.0 %>% replace_na("--")

pvalue_table <- data.frame("0.01" = p0.01, "0.1" = p0.1, "1.0" = p1.0, "2.0" = p2.0, "3.0" = p3.0, "5.0" = p5.0 )
colnames(pvalue_table) <- c(0.01, 0.1, 1.0, 2.0, 3.0, 5.0)

pfisher<- c(pfinal_0.01, pfinal_0.1, pfinal_1.0, pfinal_2.0, pfinal_3.0, pfinal_5.0)
pfisher<-scales::pvalue(pfisher)
pvalue_table[nrow(pvalue_table) + 1,]<-pfisher

row.names(pvalue_table) <- c("Block 1 - August 2nd, 2018", "Block 2 - August 6th, 2018", "Block 3 - June 14, 2019", "Block 4 - June 17, 2019", "Block 5 - June 19, 2019", "Block 6 - June 25, 2019", "Block 7 - June 27, 2019", "Block 8 - January 30, 2020", "Block 9 - February 13th, 2020", "Block 10 - September 27th, 2020", "Block 11 - October 6th, 2020", "Block 12 - November 1st, 2020", "Block 13 - June 16th, 2021", " Final P-value ")

pvalue_table %>%
  kbl(caption = "<center><strong>P-values for log-rank tests on individual doses within each experiment<center><strong>") %>%
   kable_classic(full_width = F, html_font = "Cambria", position = "center")%>%
    add_header_above(c(" "=1, "Infectious Dose (Abs600nm)" = 6)) %>%
    pack_rows("Individual Blocks", 1, 13) %>%
    pack_rows("Fisher's Combined", 14, 14)
P-values for log-rank tests on individual doses within each experiment
Infectious Dose (Abs600nm)
0.01 0.1 1 2 3 5
Individual Blocks
Block 1 - August 2nd, 2018 0.074 0.127 0.021 – – –
Block 2 - August 6th, 2018 0.018 0.001 0.027 – – –
Block 3 - June 14, 2019 0.065 <0.001 <0.001 <0.001 – –
Block 4 - June 17, 2019 0.185 0.508 <0.001 <0.001 – –
Block 5 - June 19, 2019 0.205 <0.001 0.003 <0.001 – –
Block 6 - June 25, 2019 0.006 <0.001 <0.001 <0.001 – –
Block 7 - June 27, 2019 0.650 0.008 <0.001 <0.001 – –
Block 8 - January 30, 2020 – <0.001 <0.001 <0.001 – <0.001
Block 9 - February 13th, 2020 – <0.001 <0.001 <0.001 – <0.001
Block 10 - September 27th, 2020 – – <0.001 <0.001 <0.001 <0.001
Block 11 - October 6th, 2020 – – <0.001 <0.001 <0.001 <0.001
Block 12 - November 1st, 2020 – – <0.001 <0.001 <0.001 <0.001
Block 13 - June 16th, 2021 – <0.001 <0.001 <0.001 – <0.001
Fisher’s Combined
Final P-value <0.001 <0.001 <0.001 <0.001 <0.001 <0.001
  • Conclusion: There was significant protection at all infectious doses for the secondary infection (p<0.001)

Figure 2 - E. faecalis Chronic Infection

  • Mortality over 7 days following secondary infection with P. rettgeri in flies chronically infected with E. faecalis versus control flies not given a chronic infection.

Graph

#KAPLAN MEIER
#table(Ef_Pr$primary, Ef_Pr$secondary_dose) #tells you the number of flies in each condition
#png(file="E.fae-P.rett.png", width=1200, heigh=900, res=130) #determines where the graph will be saved
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
#The graph
plot(survfit(Surv(Ef_Pr$day, Ef_Pr$status) ~ Ef_Pr$primary + Ef_Pr$secondary_dose),
     col=c("#B5E8FC", "#09B6FC", "#055df5", "#0600a8"), #line color, cycles through
     lty=c(3, 3, 3, 3, 1, 1, 1, 1), #line type
     lwd=5, #line width
     yaxt='n', #hide the y axis so I can customize it
     xaxt='n', #hide the x axis so I can customize itS
     ylab="proportion alive", #label for the y-axis
     xlab="days post-secondary infection", #label for the x-axis
     cex.lab=2.5, #font size of axis labels
     xlim=c(0,7), #determines range of the x-axis
     xaxs='i') #makes the data flush with the axis
#Y axis
axis(2, #left side
     las=2, #perpendicular to the axis
     cex.axis=1.5, #font size
     lwd=5) #line width
#X axis
axis(1, #bottom of the graph
     cex.axis=1.5, #font size
     lwd=5) #line width
box(lwd=5)
#The legend
legend("topright", #location
       bty="n", #no box around the legend
       c("PBS-P.rett 0.001", "PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.001", "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"), #labels in the legend
       col=c("#B5E8FC", "#09B6FC", "#055df5", "#0600a8"), #colors for the legend lines, cycles through
       lty=c(3, 3, 3, 3, 1, 1, 1, 1), #line type for legend lines
       lwd=5, #line width for legend lines
       cex=1) #font size

Figure 2 Panels

Figure 2A - 0.001

Ef_Pr_0.001<-subset(Ef_Pr,secondary_dose=="0.001")
#table(Ef_Pr_0.001$primary)
#THE GRAPH
#png(file="Figure2A.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Ef_Pr_0.001), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=101)", "chronic & secondary infection (n=110)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.8)

#dev.off()

Figure 2B - 0.01

Ef_Pr_0.01<-subset(Ef_Pr,secondary_dose=="0.01")
#table(Ef_Pr_0.01$primary)
#THE GRAPH
#png(file="Figure2B.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Ef_Pr_0.01), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=148)", "chronic & secondary infection (n=119)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.8)

#dev.off()

Figure 2C - 0.1

Ef_Pr_0.1<-subset(Ef_Pr,secondary_dose=="0.1")
#table(Ef_Pr_0.1$primary)
#THE GRAPH
#png(file="Figure2C.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Ef_Pr_0.1), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=160)", "chronic & secondary infection (n=115)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.8)

#dev.off()

Figure 2D - 1.0

Ef_Pr_1.0<-subset(Ef_Pr,secondary_dose=="1")
#table(Ef_Pr_1.0$primary)
#THE GRAPH
#png(file="Figure2D.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Ef_Pr_1.0), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=157)", "chronic & secondary infection (n=125)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.8)

#dev.off()

Figure 2

g <- ggdraw() +
  draw_image("Figure2A.png")

h <- ggdraw() +
  draw_image("Figure2B.png")

i <- ggdraw() +
  draw_image("Figure2C.png")

j <- ggdraw() +
  draw_image("Figure2D.png")


multiplot <- g + h + i + j
y <- multiplot + 
  plot_annotation(tag_levels = 'A')
ggsave("Fig.2.ChronicEf_Pr.pdf", y, device = "pdf")
y

Individual Blocks

Block 1

Ef_Pr_block1 <- subset(Ef_Pr, Date == "6/11/2019")

table(Ef_Pr_block1 $primary, Ef_Pr_block1 $secondary_dose)
##         
##          0.01 0.1  1
##   01_PBS   20  19 20
##   02_Ef    16  16 17
#Ef_Pr_0.001_1<-subset(Ef_Pr_block1, secondary_dose=="2")
#table(Ef_Pr_0.001_1$primary)
Ef_Pr_0.01_1<-subset(Ef_Pr_block1, secondary_dose=="0.01")
#table(Ef_Pr_0.01_1$primary)
Ef_Pr_0.1_1<-subset(Ef_Pr_block1, secondary_dose=="0.1")
#table(Ef_Pr_0.1_1$primary)
Ef_Pr_1.0_1<-subset(Ef_Pr_block1, secondary_dose=="1")
#table(Ef_Pr_1.0_1$primary)

surv_Ef_Pr_0.01_1<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.01_1)
surv_Ef_Pr_0.01_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.01_1)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 20       19     12.8      3.03      8.21
## primary=02_Ef  16       16     22.2      1.74      8.21
## 
##  Chisq= 8.2  on 1 degrees of freedom, p= 0.004
surv_Ef_Pr_0.1_1<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.1_1)
surv_Ef_Pr_0.1_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.1_1)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 19       19     14.2      1.66       5.9
## primary=02_Ef  16       15     19.8      1.18       5.9
## 
##  Chisq= 5.9  on 1 degrees of freedom, p= 0.02
surv_Ef_Pr_1.0_1<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_1.0_1)
surv_Ef_Pr_1.0_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_1.0_1)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 20       20     17.8     0.262      5.13
## primary=02_Ef  17       17     19.2     0.244      5.13
## 
##  Chisq= 5.1  on 1 degrees of freedom, p= 0.02
par(mar=c(4,5,1,1))
plot(survfit(Surv(Ef_Pr_block1$day, Ef_Pr_block1$status) ~ Ef_Pr_block1$primary + Ef_Pr_block1$secondary_dose), col=c("#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"),col=c("#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 1, 1, 1),lwd=5, cex=1)

Block 2

Ef_Pr_block2 <- subset(Ef_Pr, Date == "6/22/2018")


table(Ef_Pr_block2 $primary, Ef_Pr_block2 $secondary_dose)
##         
##          0.001 0.01 0.1  1
##   01_PBS    28   18  29 30
##   02_Ef     33   18  19 29
Ef_Pr_0.001_2<-subset(Ef_Pr_block2, secondary_dose=="0.001")
table(Ef_Pr_0.001_2$primary)
## 
## 01_PBS  02_Ef 
##     28     33
Ef_Pr_0.01_2<-subset(Ef_Pr_block2, secondary_dose=="0.01")
#table(Ef_Pr_0.01_2$primary)
Ef_Pr_0.1_2<-subset(Ef_Pr_block2, secondary_dose=="0.1")
#table(Ef_Pr_0.1_2$primary)
Ef_Pr_1.0_2<-subset(Ef_Pr_block2, secondary_dose=="1")
#table(Ef_Pr_1.0_2$primary)

surv_Ef_Pr_0.001_2<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.001_2)
surv_Ef_Pr_0.001_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.001_2)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 28       13     10.9     0.406      1.01
## primary=02_Ef  33       11     13.1     0.337      1.01
## 
##  Chisq= 1  on 1 degrees of freedom, p= 0.3
surv_Ef_Pr_0.01_2<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.01_2)
surv_Ef_Pr_0.01_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.01_2)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 18       13     8.44      2.46      6.03
## primary=02_Ef  18        6    10.56      1.97      6.03
## 
##  Chisq= 6  on 1 degrees of freedom, p= 0.01
surv_Ef_Pr_0.1_2<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.1_2)
surv_Ef_Pr_0.1_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.1_2)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 29       23     17.1      2.01      7.88
## primary=02_Ef  19        8     13.9      2.48      7.88
## 
##  Chisq= 7.9  on 1 degrees of freedom, p= 0.005
surv_Ef_Pr_1.0_2<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_1.0_2)
surv_Ef_Pr_1.0_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_1.0_2)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 30       30     27.5     0.235      5.56
## primary=02_Ef  29       26     28.5     0.226      5.56
## 
##  Chisq= 5.6  on 1 degrees of freedom, p= 0.02
par(mar=c(4,5,1,1))
plot(survfit(Surv(Ef_Pr_block2$day, Ef_Pr_block2$status) ~ Ef_Pr_block2$primary + Ef_Pr_block2$secondary_dose), col=c("#B5E8FC",  "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.001", "PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.001", "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)

Block 3

Ef_Pr_block3 <- subset(Ef_Pr, Date == "6/25/2018")

table(Ef_Pr_block3 $primary, Ef_Pr_block3 $secondary_dose)
##         
##          0.001 0.01 0.1  1
##   01_PBS    25   27  29 37
##   02_Ef     23   15  25 19
Ef_Pr_0.001_3<-subset(Ef_Pr_block3, secondary_dose=="0.001")
table(Ef_Pr_0.001_3$primary)
## 
## 01_PBS  02_Ef 
##     25     23
Ef_Pr_0.01_3<-subset(Ef_Pr_block3, secondary_dose=="0.01")
#table(Ef_Pr_0.01_3$primary)
Ef_Pr_0.1_3<-subset(Ef_Pr_block3, secondary_dose=="0.1")
#table(Ef_Pr_0.1_3$primary)
Ef_Pr_1.0_3<-subset(Ef_Pr_block3, secondary_dose=="1")
#table(Ef_Pr_1.0_3$primary)

surv_Ef_Pr_0.001_3<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.001_3)
surv_Ef_Pr_0.001_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.001_3)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 25       18     14.4     0.912      3.13
## primary=02_Ef  23       11     14.6     0.897      3.13
## 
##  Chisq= 3.1  on 1 degrees of freedom, p= 0.08
surv_Ef_Pr_0.01_3<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.01_3)
surv_Ef_Pr_0.01_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.01_3)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 27       21     15.8      1.68      7.31
## primary=02_Ef  15        5     10.2      2.62      7.31
## 
##  Chisq= 7.3  on 1 degrees of freedom, p= 0.007
surv_Ef_Pr_0.1_3<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.1_3)
surv_Ef_Pr_0.1_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.1_3)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 29       21     20.6   0.00806    0.0384
## primary=02_Ef  25       18     18.4   0.00902    0.0384
## 
##  Chisq= 0  on 1 degrees of freedom, p= 0.8
surv_Ef_Pr_1.0_3<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_1.0_3)
surv_Ef_Pr_1.0_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_1.0_3)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 37       36     35.7   0.00290     0.235
## primary=02_Ef  19       18     18.3   0.00564     0.235
## 
##  Chisq= 0.2  on 1 degrees of freedom, p= 0.6
par(mar=c(4,5,1,1))
plot(survfit(Surv(Ef_Pr_block3$day, Ef_Pr_block3$status) ~ Ef_Pr_block3$primary + Ef_Pr_block3$secondary_dose), col=c("#B5E8FC",  "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.001", "PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.001", "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)

Block 4

Not sure why the data stops after 3 days, Frank did this one

Ef_Pr_block4 <- subset(Ef_Pr, Date == "7/19/2018")

table(Ef_Pr_block4 $primary, Ef_Pr_block4 $secondary_dose)
##         
##          0.001 0.01 0.1  1
##   01_PBS     9   18  18 10
##   02_Ef     23   27  16 17
Ef_Pr_0.001_4<-subset(Ef_Pr_block4, secondary_dose=="0.001")
#table(Ef_Pr_0.001_4$primary)
Ef_Pr_0.01_4<-subset(Ef_Pr_block4, secondary_dose=="0.01")
#table(Ef_Pr_0.01_4$primary)
Ef_Pr_0.1_4<-subset(Ef_Pr_block4, secondary_dose=="0.1")
#table(Ef_Pr_0.1_4$primary)
Ef_Pr_1.0_4<-subset(Ef_Pr_block4, secondary_dose=="1")
#table(Ef_Pr_1.0_4$primary)

surv_Ef_Pr_0.001_4<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.001_4)
surv_Ef_Pr_0.001_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.001_4)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS  9        3     2.56    0.0743     0.123
## primary=02_Ef  23        6     6.44    0.0296     0.123
## 
##  Chisq= 0.1  on 1 degrees of freedom, p= 0.7
surv_Ef_Pr_0.01_4<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.01_4)
surv_Ef_Pr_0.01_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.01_4)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 18       12       12  7.93e-05   0.00026
## primary=02_Ef  27       16       16  5.98e-05   0.00026
## 
##  Chisq= 0  on 1 degrees of freedom, p= 1
surv_Ef_Pr_0.1_4<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.1_4)
surv_Ef_Pr_0.1_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.1_4)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 18        9     8.31    0.0574     0.147
## primary=02_Ef  16        9     9.69    0.0492     0.147
## 
##  Chisq= 0.1  on 1 degrees of freedom, p= 0.7
surv_Ef_Pr_1.0_4<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_1.0_4)
surv_Ef_Pr_1.0_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_1.0_4)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 10        9     8.79   0.00516    0.0292
## primary=02_Ef  17       17    17.21   0.00263    0.0292
## 
##  Chisq= 0  on 1 degrees of freedom, p= 0.9
par(mar=c(4,5,1,1))
plot(survfit(Surv(Ef_Pr_block4$day, Ef_Pr_block4$status) ~ Ef_Pr_block4$primary + Ef_Pr_block4$secondary_dose), col=c("#B5E8FC",  "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.001", "PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.001", "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)

Block 5

Ef_Pr_block5 <- subset(Ef_Pr, Date == "6/3/2019")

table(Ef_Pr_block5 $primary, Ef_Pr_block5 $secondary_dose)
##         
##          0.001 0.01 0.1  1
##   01_PBS    18   21  19 17
##   02_Ef     12   13   6 12
Ef_Pr_0.001_5<-subset(Ef_Pr_block5, secondary_dose=="0.001")
table(Ef_Pr_0.001_5$primary)
## 
## 01_PBS  02_Ef 
##     18     12
Ef_Pr_0.01_5<-subset(Ef_Pr_block5, secondary_dose=="0.01")
#table(Ef_Pr_0.01_5$primary)
Ef_Pr_0.1_5<-subset(Ef_Pr_block5, secondary_dose=="0.1")
#table(Ef_Pr_0.1_5$primary)
Ef_Pr_1.0_5<-subset(Ef_Pr_block5, secondary_dose=="1")
#table(Ef_Pr_1.0_5$primary)

surv_Ef_Pr_0.001_5<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.001_5)
surv_Ef_Pr_0.001_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.001_5)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 18       13    10.65     0.518      1.39
## primary=02_Ef  12        6     8.35     0.661      1.39
## 
##  Chisq= 1.4  on 1 degrees of freedom, p= 0.2
surv_Ef_Pr_0.01_5<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.01_5)
surv_Ef_Pr_0.01_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.01_5)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 21       16     12.9     0.759      2.31
## primary=02_Ef  13        7     10.1     0.965      2.31
## 
##  Chisq= 2.3  on 1 degrees of freedom, p= 0.1
surv_Ef_Pr_0.1_5<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.1_5)
surv_Ef_Pr_0.1_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.1_5)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 19       15    14.05    0.0636     0.502
## primary=02_Ef   6        4     4.95    0.1806     0.502
## 
##  Chisq= 0.5  on 1 degrees of freedom, p= 0.5
surv_Ef_Pr_1.0_5<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_1.0_5)
surv_Ef_Pr_1.0_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_1.0_5)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 17       17     14.7     0.375      6.35
## primary=02_Ef  12       12     14.3     0.383      6.35
## 
##  Chisq= 6.3  on 1 degrees of freedom, p= 0.01
par(mar=c(4,5,1,1))
plot(survfit(Surv(Ef_Pr_block5$day, Ef_Pr_block5$status) ~ Ef_Pr_block5$primary + Ef_Pr_block5$secondary_dose), col=c("#B5E8FC",  "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.001", "PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.001", "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)

Block 6

Ef_Pr_block6 <- subset(Ef_Pr, Date == "6/5/2019")

table(Ef_Pr_block6 $primary, Ef_Pr_block6 $secondary_dose)
##         
##          0.001 0.01 0.1  1
##   01_PBS    21   21  22 20
##   02_Ef     19   16  16 17
Ef_Pr_0.001_6<-subset(Ef_Pr_block6, secondary_dose=="0.001")
table(Ef_Pr_0.001_6$primary)
## 
## 01_PBS  02_Ef 
##     21     19
Ef_Pr_0.01_6<-subset(Ef_Pr_block6, secondary_dose=="0.01")
#table(Ef_Pr_0.01_6$primary)
Ef_Pr_0.1_6<-subset(Ef_Pr_block6, secondary_dose=="0.1")
#table(Ef_Pr_0.1_6$primary)
Ef_Pr_1.0_6<-subset(Ef_Pr_block6, secondary_dose=="1")
#table(Ef_Pr_1.0_6$primary)

surv_Ef_Pr_0.001_6<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.001_6)
surv_Ef_Pr_0.001_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.001_6)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 21       13     8.18      2.84      6.64
## primary=02_Ef  19        4     8.82      2.63      6.64
## 
##  Chisq= 6.6  on 1 degrees of freedom, p= 0.01
surv_Ef_Pr_0.01_6<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.01_6)
surv_Ef_Pr_0.01_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.01_6)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 21       10     7.63     0.738      1.76
## primary=02_Ef  16        4     6.37     0.883      1.76
## 
##  Chisq= 1.8  on 1 degrees of freedom, p= 0.2
surv_Ef_Pr_0.1_6<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.1_6)
surv_Ef_Pr_0.1_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.1_6)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 22       10    10.96    0.0834     0.246
## primary=02_Ef  16        9     8.04    0.1137     0.246
## 
##  Chisq= 0.2  on 1 degrees of freedom, p= 0.6
surv_Ef_Pr_1.0_6<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_1.0_6)
surv_Ef_Pr_1.0_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_1.0_6)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 20       17     16.7   0.00591     0.026
## primary=02_Ef  17       16     16.3   0.00604     0.026
## 
##  Chisq= 0  on 1 degrees of freedom, p= 0.9
par(mar=c(4,5,1,1))
plot(survfit(Surv(Ef_Pr_block6$day, Ef_Pr_block6$status) ~ Ef_Pr_block6$primary + Ef_Pr_block6$secondary_dose), col=c("#B5E8FC",  "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.001", "PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.001", "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)

Block 7

Ef_Pr_block7 <- subset(Ef_Pr, Date == "7/1/2019")

table(Ef_Pr_block7 $primary, Ef_Pr_block7 $secondary_dose)
##         
##          0.01 0.1  1
##   01_PBS   23  24 23
##   02_Ef    14  17 14
#Ef_Pr_0.001_7<-subset(Ef_Pr_block3, secondary_dose=="0.001")
#table(Ef_Pr_0.001_7$primary)
Ef_Pr_0.01_7<-subset(Ef_Pr_block3, secondary_dose=="0.01")
#table(Ef_Pr_0.01_7$primary)
Ef_Pr_0.1_7<-subset(Ef_Pr_block3, secondary_dose=="0.1")
#table(Ef_Pr_0.1_7$primary)
Ef_Pr_1.0_7<-subset(Ef_Pr_block3, secondary_dose=="1")
#table(Ef_Pr_1.0_7$primary)

#surv_Ef_Pr_0.001_7<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.001_7)
#surv_Ef_Pr_0.001_7

surv_Ef_Pr_0.01_7<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.01_7)
surv_Ef_Pr_0.01_7
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.01_7)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 27       21     15.8      1.68      7.31
## primary=02_Ef  15        5     10.2      2.62      7.31
## 
##  Chisq= 7.3  on 1 degrees of freedom, p= 0.007
surv_Ef_Pr_0.1_7<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.1_7)
surv_Ef_Pr_0.1_7
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.1_7)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 29       21     20.6   0.00806    0.0384
## primary=02_Ef  25       18     18.4   0.00902    0.0384
## 
##  Chisq= 0  on 1 degrees of freedom, p= 0.8
surv_Ef_Pr_1.0_7<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_1.0_7)
surv_Ef_Pr_1.0_7
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_1.0_7)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 37       36     35.7   0.00290     0.235
## primary=02_Ef  19       18     18.3   0.00564     0.235
## 
##  Chisq= 0.2  on 1 degrees of freedom, p= 0.6
par(mar=c(4,5,1,1))
plot(survfit(Surv(Ef_Pr_block7$day, Ef_Pr_block7$status) ~ Ef_Pr_block7$primary + Ef_Pr_block7$secondary_dose), col=c("#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"),col=c("#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3,  1, 1, 1),lwd=5, cex=1)

Statistical Tests

  • In order test whether there was protection provided by chronic infection for each secondary infectious dose, p-values from the log rank test for each date was combined using the Fisher’s test (experiments on all dates showed higher survival in conditions with chronic infection). Final p-value is reported in the text.
# Grab p-values from each individual logrank test (results shown above in each block)
pvalue_Ef_Pr_0.01_1<-pchisq(surv_Ef_Pr_0.01_1$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.1_1<-pchisq(surv_Ef_Pr_0.1_1$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_1.0_1<-pchisq(surv_Ef_Pr_1.0_1$chisq, df=1, lower.tail=F)

pvalue_Ef_Pr_0.001_2<-pchisq(surv_Ef_Pr_0.001_2$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.01_2<-pchisq(surv_Ef_Pr_0.01_2$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.1_2<-pchisq(surv_Ef_Pr_0.1_2$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_1.0_2<-pchisq(surv_Ef_Pr_1.0_2$chisq, df=1, lower.tail=F)

pvalue_Ef_Pr_0.001_3<-pchisq(surv_Ef_Pr_0.001_3$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.01_3<-pchisq(surv_Ef_Pr_0.01_3$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.1_3<-pchisq(surv_Ef_Pr_0.1_3$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_1.0_3<-pchisq(surv_Ef_Pr_1.0_3$chisq, df=1, lower.tail=F)

pvalue_Ef_Pr_0.001_4<-pchisq(surv_Ef_Pr_0.001_4$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.01_4<-pchisq(surv_Ef_Pr_0.01_4$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.1_4<-pchisq(surv_Ef_Pr_0.1_4$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_1.0_4<-pchisq(surv_Ef_Pr_1.0_4$chisq, df=1, lower.tail=F)

pvalue_Ef_Pr_0.001_5<-pchisq(surv_Ef_Pr_0.001_5$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.01_5<-pchisq(surv_Ef_Pr_0.01_5$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.1_5<-pchisq(surv_Ef_Pr_0.1_5$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_1.0_5<-pchisq(surv_Ef_Pr_1.0_5$chisq, df=1, lower.tail=F)

pvalue_Ef_Pr_0.001_6<-pchisq(surv_Ef_Pr_0.001_6$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.01_6<-pchisq(surv_Ef_Pr_0.01_6$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.1_6<-pchisq(surv_Ef_Pr_0.1_6$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_1.0_6<-pchisq(surv_Ef_Pr_1.0_6$chisq, df=1, lower.tail=F)

pvalue_Ef_Pr_0.01_7<-pchisq(surv_Ef_Pr_0.01_7$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.1_7<-pchisq(surv_Ef_Pr_0.1_7$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_1.0_7<-pchisq(surv_Ef_Pr_1.0_7$chisq, df=1, lower.tail=F)
# Use all p-values from a given condition to get Fisher's combined test p-value# Use all p-values from a given condition to get Fisher's combined test p-value
pvalues_Ef_Pr_0.001 <- c(pvalue_Ef_Pr_0.001_2,pvalue_Ef_Pr_0.001_3,pvalue_Ef_Pr_0.001_4,pvalue_Ef_Pr_0.001_5,pvalue_Ef_Pr_0.001_6)
test.stat_Ef_Pr_0.001 <- -2*sum(log(pvalues_Ef_Pr_0.001))
pfinal_Ef_Pr_0.001<-pchisq(test.stat_Ef_Pr_0.001, df=2*length(pvalues_Ef_Pr_0.001), lower.tail=F)

pvalues_Ef_Pr_0.01 <- c(pvalue_Ef_Pr_0.01_1,pvalue_Ef_Pr_0.01_2,pvalue_Ef_Pr_0.01_3,pvalue_Ef_Pr_0.01_4,pvalue_Ef_Pr_0.01_5,pvalue_Ef_Pr_0.01_6, pvalue_Ef_Pr_0.01_7)
test.stat_Ef_Pr_0.01 <- -2*sum(log(pvalues_Ef_Pr_0.01))
pfinal_Ef_Pr_0.01<-pchisq(test.stat_Ef_Pr_0.01, df=2*length(pvalues_Ef_Pr_0.01), lower.tail=F)


pvalues_Ef_Pr_0.1 <- c(pvalue_Ef_Pr_0.1_1,pvalue_Ef_Pr_0.1_2,pvalue_Ef_Pr_0.1_3,pvalue_Ef_Pr_0.1_4,pvalue_Ef_Pr_0.1_5,pvalue_Ef_Pr_0.1_6,pvalue_Ef_Pr_0.1_7)
test.stat_Ef_Pr_0.1 <- -2*sum(log(pvalues_Ef_Pr_0.1))
pfinal_Ef_Pr_0.1<-pchisq(test.stat_Ef_Pr_0.1, df=2*length(pvalues_Ef_Pr_0.1), lower.tail=F)

pvalues_Ef_Pr_1.0 <- c(pvalue_Ef_Pr_1.0_1,pvalue_Ef_Pr_1.0_2,pvalue_Ef_Pr_1.0_3,pvalue_Ef_Pr_1.0_4,pvalue_Ef_Pr_1.0_5,pvalue_Ef_Pr_1.0_6,pvalue_Ef_Pr_1.0_7)
test.stat_Ef_Pr_1.0 <- -2*sum(log(pvalues_Ef_Pr_1.0))
pfinal_Ef_Pr_1.0<-pchisq(test.stat_Ef_Pr_1.0, df=2*length(pvalues_Ef_Pr_1.0), lower.tail=F)
####  Generate Table of final results
pvalues_Ef_Pr_0.001_table <- c(NA,pvalue_Ef_Pr_0.001_2,pvalue_Ef_Pr_0.001_3,pvalue_Ef_Pr_0.001_4,pvalue_Ef_Pr_0.001_5,pvalue_Ef_Pr_0.001_6,NA)
p_Ef_Pr_0.001<-scales::pvalue(pvalues_Ef_Pr_0.001_table)
p_Ef_Pr_0.001 <- p_Ef_Pr_0.001 %>% replace_na("--")

pvalues_Ef_Pr_0.01_table <- c(pvalue_Ef_Pr_0.01_1,pvalue_Ef_Pr_0.01_2,pvalue_Ef_Pr_0.01_3,pvalue_Ef_Pr_0.01_4,pvalue_Ef_Pr_0.01_5,pvalue_Ef_Pr_0.01_6,pvalue_Ef_Pr_0.01_7)
p_Ef_Pr_0.01<-scales::pvalue(pvalues_Ef_Pr_0.01_table)

pvalues_Ef_Pr_0.1_table <- c(pvalue_Ef_Pr_0.1_1,pvalue_Ef_Pr_0.1_2,pvalue_Ef_Pr_0.1_3,pvalue_Ef_Pr_0.1_4,pvalue_Ef_Pr_0.1_5,pvalue_Ef_Pr_0.1_6,pvalue_Ef_Pr_0.1_7)
p_Ef_Pr_0.1<-scales::pvalue(pvalues_Ef_Pr_0.1_table)

pvalues_Ef_Pr_1.0_table <- c(pvalue_Ef_Pr_1.0_1,pvalue_Ef_Pr_1.0_2,pvalue_Ef_Pr_1.0_3,pvalue_Ef_Pr_1.0_4,pvalue_Ef_Pr_1.0_5,pvalue_Ef_Pr_1.0_6,pvalue_Ef_Pr_1.0_7)
p_Ef_Pr_1.0<-scales::pvalue(pvalues_Ef_Pr_1.0_table)



pvalue_table_Ef_Pr <- data.frame("0.001" = p_Ef_Pr_0.001, "0.01" = p_Ef_Pr_0.01, "0.1" = p_Ef_Pr_0.1, "1.0" = p_Ef_Pr_1.0)
colnames(pvalue_table_Ef_Pr) <- c(0.001,0.01, 0.1, 1.0)

pfisher_Ef_Pr<- c(pfinal_Ef_Pr_0.001, pfinal_Ef_Pr_0.01, pfinal_Ef_Pr_0.1, pfinal_Ef_Pr_1.0)
pfisher_Ef_Pr<-scales::pvalue(pfisher_Ef_Pr)
pvalue_table_Ef_Pr[nrow(pvalue_table_Ef_Pr) + 1,]<-pfisher_Ef_Pr

row.names(pvalue_table_Ef_Pr) <- c("Block 1 - June 11th, 2018", "Block 2 - June 22nd, 2018", "Block 3 - June 25th, 2018", "Block 4 - July 19th, 2018", "Block 5 - June 3, 2019", "Block 6 - June 5, 2019", "Block 7 - July 1st, 2019", " Final P-value ")

pvalue_table_Ef_Pr %>%
  kbl(caption = "<center><strong>P-values for log-rank tests on individual doses within each experiment<center><strong>") %>%
   kable_classic(full_width = F, html_font = "Cambria", position = "center")%>%
    add_header_above(c(" "=1, "Infectious Dose (Abs600nm)" = 4)) %>%
    pack_rows("Individual Blocks", 1, 7) %>%
    pack_rows("Fisher's Combined", 8, 8)
P-values for log-rank tests on individual doses within each experiment
Infectious Dose (Abs600nm)
0.001 0.01 0.1 1
Individual Blocks
Block 1 - June 11th, 2018 – 0.004 0.015 0.023
Block 2 - June 22nd, 2018 0.315 0.014 0.005 0.018
Block 3 - June 25th, 2018 0.077 0.007 0.845 0.628
Block 4 - July 19th, 2018 0.725 0.987 0.702 0.864
Block 5 - June 3, 2019 0.239 0.129 0.479 0.012
Block 6 - June 5, 2019 0.010 0.184 0.620 0.872
Block 7 - July 1st, 2019 – 0.007 0.845 0.628
Fisher’s Combined
Final P-value 0.028 <0.001 0.064 0.020
  • Conclusion: Chronic infection with E. faecalis offered significat protection against three out of the four doses. However, the one dose that did not show significant protection (Abs600nm = 0.1) had a borderline p-value = 0.06

Q2 Does having a primary infection influence resistance to secondary infection with P. rettgeri?

Study design:

  • Bacterial load of P. rettgeri was determined at 10 hours post-secondary infection
  • For S. marcescens-P. rettgeri, the load was determined by isolating pooled DNA from 5 flies for each condition and using qPCR primer specific for a P. rettgeri gene to assess the relative number of bacteria. This was converted to approximate bacterial load by using a conversion factor based on qPCR and CFU results for flies injected with 20nL of 0.1 Abs600nm of P. rettgeri and then immediately homogenized.
  • For E. faecalis - P. rettgeri, the P. rettgeri load was determined using a Colony Forming Unit (CFU) assay using antibiotics to select for P. rettgeri and kill E. faecalis.
  • All analyses are with linear modeling followed by anova, with a shapiro test on residuals to assess normality.
  • Conclusions
    • As expected based on previously published results (Chambers et al 2019), infectious dose of the secondary P. rettgeri infection significantly impacted bacterial load 10 hours after infection (p<0.0001)
    • Chronic infection with S. marcescens and E. faecalis both resulted in lowered calculated bacterial load post secondary infection (p<0.0001).
    • There was no significant interaction between chronic infection and secondary infectious dose for either S. marcescens or E. faecalis chronic infection (p=0.30 and 0.78 respectively)
    • Residuals pass normality test (Shapiro-Wilk Normality test: p=0.72 and 0.79 respectively)

S. marcescens chronic infection

  • The calculated bacterial load after secondary infection with P. rettgeri for flies chronically infected with S. marcescens versus control flies.
  • As expected based on previously published results (Chambers et al 2019), infectious dose of the secondary P. rettgeri infection significantly impacted bacterial load 10 hours after infection (secondary_dose, p<0.0001)
  • Primary infection with S. marcescens resulted in lowered calculated bacterial load post secondary infection (primary, p<0.0001).
  • There was no significant interaction between primary infection and secondary infectious dose (p=0.30)
  • Residuals pass normality test (Shapiro-Wilk Normality test: p=0.72)

Graph

#png(file="Figure3A.png", width=3200, height=1800, res=300)
ggplot(Sm_bacterial_load, aes(x = as.character(secondary_dose), y = log(load_10hr), fill = primary)) +
  scale_fill_manual(values=c("white", "#3399FF"), 
                       name=" ",
                       breaks=c("01_PBS","02_Sm"),
                       labels=c("control", "chronically infected")) +
  geom_boxplot() +
  #ggtitle("S. marcescens chronic infection") +
  theme_classic(base_size = 24) +
  xlab("secondary dose (Abs600nm)") +
  ylab("log(calculated bacterial load)") +
  guides(fill=guide_legend(title=NULL)) +
  geom_point(color = "black", size = 1.3, position=position_jitterdodge(jitter.width=0.2))

#dev.off()

Statistics

S.m_resistance <- lm(log(load_10hr) ~ rep + secondary_dose + primary + secondary_dose:primary, data = Sm_bacterial_load)
S.m_resistance_anova <- anova(S.m_resistance)
S.m_resistance_anova
## Analysis of Variance Table
## 
## Response: log(load_10hr)
##                        Df  Sum Sq Mean Sq F value    Pr(>F)    
## rep                     8  56.729   7.091  2.8793   0.01059 *  
## secondary_dose          1 145.561 145.561 59.1044 6.543e-10 ***
## primary                 1  94.077  94.077 38.1995 1.328e-07 ***
## secondary_dose:primary  1   2.651   2.651  1.0765   0.30467    
## Residuals              48 118.213   2.463                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
shapiro.test(resid(S.m_resistance))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(S.m_resistance)
## W = 0.98595, p-value = 0.7192

E. faecalis chronic infection

  • The bacterial load after secondary infection with P. rettgeri for flies chronically infected with E. faecalis versus control flies.
  • As expected based on previously published results (Chambers et al 2019), infectious dose of the secondary P. rettgeri infection significantly impacted bacterial load 10 hours after infection (secondary_dose, p<0.0001)
  • Primary infection with E. faecalis resulted in lowered CFU/fly post secondary infection (primary, p = 0.0001).
  • There was no significant interaction between primary infection and secondary infectious dose (p=0.78)
  • Residuals pass normality test (Shapiro-Wilk Normality test: p=0.79)

Graph

#BACTERIAL LOAD
Ef_bacterial_load_10hr <- subset(Ef_bacterial_load, time == 10)

#png(file="Figure3B.png",width=3200, height=1800, res=300)
ggplot(Ef_bacterial_load_10hr, aes(x = as.character(secondary_dose), y = log(load), fill = primary)) +
  scale_fill_manual(values=c("white", "#3399FF"), 
                       name=" ",
                       breaks=c("01_PBS","02_Efae"),
                       labels=c("control", "chronically infected")) +
  geom_boxplot() +
  #ggtitle("E. faecalis chronic infection") +
  theme_classic(base_size = 24) +
  xlab("secondary dose (Abs600nm)") +
  ylab("log(CFU/fly)") +
  guides(fill=guide_legend(title=NULL)) +
  geom_point(color = "black", size = 1.3, position=position_jitterdodge(jitter.width=0.2))

#dev.off()

Statistics

E.fae_resistance <- lm(log(load) ~ replicate + secondary_dose + primary + secondary_dose:primary, data = Ef_bacterial_load_10hr)
E.fae_resistance_anova <- anova(E.fae_resistance)
E.fae_resistance_anova
## Analysis of Variance Table
## 
## Response: log(load)
##                         Df Sum Sq Mean Sq F value    Pr(>F)    
## replicate                5 231.74   46.35  15.403 3.504e-13 ***
## secondary_dose           1 799.17  799.17 265.589 < 2.2e-16 ***
## primary                  1  35.09   35.09  11.660 0.0007463 ***
## secondary_dose:primary   1   0.23    0.23   0.075 0.7843830    
## Residuals              247 743.23    3.01                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
shapiro.test(resid(E.fae_resistance))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(E.fae_resistance)
## W = 0.99623, p-value = 0.7984

Figure 3

k <- ggdraw() +
  ggtitle('P. rettgeri') +
  draw_image("Figure3A.png")

l <- ggdraw() +
  draw_image("Figure3B.png")

m <- ggdraw() +
  draw_image("Figure3C.png")

n <- ggdraw() +
  draw_image("Figure3D.png")

multiplot3 <- k + m + l + n

x<-multiplot3 
 
  #plot_annotation(tag_levels = 'A')
 
ggsave("Fig.3.ResistanceTolerance.pdf", x, device = "pdf")
x

Q3 Does the protection conferred by primary infection with S. marcescens protect against infections that cause very high mortality?

Study design:

  • Since we saw such robust protection with S. marcescens, we wanted to determine whether it also protected against a more lethal infection with P. sneebia.
  • Initially analyzed using Cox Proportional Hazards mixed effects model with date of injection, primary injection and secondary dose as fixed effects and random effect terms accounting for experimental structure (e.g. date of injection, housing of groups of flies in vials). However, this model violated the assumptions of Cox Proportional Hazards. + Therefore we switched to using log-rank pairwise comparisons to determine whether chronic infection significantly impacted survival of secondary infection at an individual dose.
    • LogRank comparisons were done on data from each date individually
    • If experiments on all dates showed higher survival in conditions with chronic infection, p-values were combined using Fisher’s test and only the final p-value is reported in the text.
    • If there was a variable effect of chronic infection, the number of individual experiments with significant effects is reported.
      • Conclusions
    • S. marcescens chronic infection provided significant protection during secondary infection with P. sneebia at all infectious doses (p<0.001)

Total Graph

#table(Sm_Ps$primary, Sm_Ps$secondary_dose)

#png(file="Sm_Ps.png",  width = 1700, height = 1500, res = 300)
#par(mar=c(4,4,0.1,0.1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(Sm_Ps$hour, Sm_Ps$status) ~ Sm_Ps$primary + Sm_Ps$secondary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=2.6, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)

legend("bottomleft", bty="n" , c("PBS-PBS", "PBS-0.001", "PBS-0.01", "PBS-0.1", "PBS-1.0", "S.m-PBS", "S.m-0.001", "S.m-0.01", "S.m-0.1", "S.m-1.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=4, cex=1)

Appendix Figure 1B - Controls

Sm_Ps_PBS<-subset(Sm_Ps,secondary_dose=="0")
#table(Sm_Ps_PBS$primary)

#THE GRAPH
#png(file="AFigure1B.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Ps_PBS), col=c("gray"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("bottomleft",inset=0, bty="n" , c("double injection control (n=59)", "chronic infection only (n=59)"),col=c("gray"),lty=c(3, 1),lwd=5, cex=1)

#dev.off()

Individual Figure Panels

Figure 4A - 0.001

Sm_Ps_0.001<-subset(Sm_Ps,secondary_dose=="0.001")
#table(Sm_Ps_0.001$primary)
#THE GRAPH
#png(file="Figure4A.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(hour, status) ~ primary + secondary_dose, data=Sm_Ps_0.001), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="hours post-secondary infection", cex.lab=2, xlim=c(0,72), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("bottomleft",inset=0, bty="n" , c("secondary infection only (n=82)", "chronic & secondary infection (n=66)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.6)

#dev.off()

Figure 4B - 0.01

Sm_Ps_0.01<-subset(Sm_Ps,secondary_dose=="0.01")
#table(Sm_Ps_0.01$primary)
#THE GRAPH
#png(file="Figure4B.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(hour, status) ~ primary + secondary_dose, data=Sm_Ps_0.01), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="hours post-secondary infection", cex.lab=2, xlim=c(0,72), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("bottomleft",inset=0, bty="n" , c("secondary infection only (n=87)", "chronic & secondary infection (n=80)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.6)

#dev.off()

Figure 4C - 0.1

Sm_Ps_0.1<-subset(Sm_Ps,secondary_dose=="0.1")
#table(Sm_Ps_0.1$primary)
#THE GRAPH
#png(file="Figure4C.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(hour, status) ~ primary + secondary_dose, data=Sm_Ps_0.1), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="hours post-secondary infection", cex.lab=2, xlim=c(0,72), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=80)", "chronic & secondary infection (n=81)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.6)

#dev.off()

Figure 4D - 1.0

Sm_Ps_1.0<-subset(Sm_Ps,secondary_dose=="1")
#table(Sm_Ps_1.0$primary)
#png(file="Figure4D.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(hour, status) ~ primary + secondary_dose, data=Sm_Ps_1.0), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="hours post-secondary infection", cex.lab=2, xlim=c(0,72), xaxs='i') 
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=75)", "chronic & secondary infection (n=81)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.6)

#dev.off()

Figure 4

m <- ggdraw() +
  draw_image("Figure4A.png")

n <- ggdraw() +
  draw_image("Figure4B.png")

o <- ggdraw() +
  ggtitle('P. rettgeri') +
  draw_image("Figure4C.png")

p <- ggdraw() +
  draw_image("Figure4D.png")


multiplot3 <- m + n + o + p

w<-multiplot3 +
 
  plot_annotation(tag_levels = 'A')
 
ggsave("Fig.4.Psneebia.pdf", w, device = "pdf")
w

Individual Blocks

Block 1

Sm_Ps_block1 <- subset(Sm_Ps, date == "10/1/2020")
table(Sm_Ps_block1$primary, Sm_Ps_block1$secondary_dose)
##      
##        0 0.001 0.01 0.1  1
##   PBS 19    30   29  20 30
##   S.m 19    16   24  29 25
#subset by secondary dose
Sm_Ps_0.001_1<-subset(Sm_Ps_block1, secondary_dose=="0.001")
Sm_Ps_0.01_1<-subset(Sm_Ps_block1, secondary_dose=="0.01")
Sm_Ps_0.1_1<-subset(Sm_Ps_block1, secondary_dose=="0.1")
Sm_Ps_1.0_1<-subset(Sm_Ps_block1, secondary_dose=="1")

surv_Sm_Ps_0.001_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.001_1)
surv_Sm_Ps_0.001_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.001_1)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30       21     18.3     0.388      1.13
## primary=S.m 16       10     12.7     0.562      1.13
## 
##  Chisq= 1.1  on 1 degrees of freedom, p= 0.3
surv_Sm_Ps_0.01_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.01_1)
surv_Sm_Ps_0.01_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.01_1)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 29       29     14.7     13.87      29.7
## primary=S.m 24       23     37.3      5.47      29.7
## 
##  Chisq= 29.7  on 1 degrees of freedom, p= 5e-08
surv_Sm_Ps_0.1_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.1_1)
surv_Sm_Ps_0.1_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.1_1)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 20       20      6.8     25.65      45.2
## primary=S.m 29       29     42.2      4.13      45.2
## 
##  Chisq= 45.2  on 1 degrees of freedom, p= 2e-11
surv_Sm_Ps_1.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_1.0_1)
surv_Sm_Ps_1.0_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_1.0_1)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30       30     18.3      7.51      19.9
## primary=S.m 25       25     36.7      3.74      19.9
## 
##  Chisq= 19.9  on 1 degrees of freedom, p= 8e-06
plot(survfit(Surv(Sm_Ps_block1$hour, Sm_Ps_block1$status) ~ Sm_Ps_block1$primary + Sm_Ps_block1$secondary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=2.3, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)

legend("bottomleft", bty="n" , c("PBS-PBS", "0.001", "PBS-0.01", "PBS-0.1", "PBS-1.0", "S.m-PBS", "S.m-0.001", "S.m-0.01", "S.m-0.1", "S.m-1.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=4, cex=1)

Block 2

Sm_Ps_block2 <- subset(Sm_Ps, date == "10/10/2020")
table(Sm_Ps_block2$primary, Sm_Ps_block2$secondary_dose)
##      
##        0 0.001 0.01 0.1  1
##   PBS 20    30   30  30 30
##   S.m 20    30   30  30 30
#subset by secondary dose
Sm_Ps_0.001_2<-subset(Sm_Ps_block2, secondary_dose=="0.001")
Sm_Ps_0.01_2<-subset(Sm_Ps_block2, secondary_dose=="0.01")
Sm_Ps_0.1_2<-subset(Sm_Ps_block2, secondary_dose=="0.1")
Sm_Ps_1.0_2<-subset(Sm_Ps_block2, secondary_dose=="1")

surv_Sm_Ps_0.001_2<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.001_2)
surv_Sm_Ps_0.001_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.001_2)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30       25       16      5.07      10.5
## primary=S.m 30       13       22      3.69      10.5
## 
##  Chisq= 10.5  on 1 degrees of freedom, p= 0.001
surv_Sm_Ps_0.01_2<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.01_2)
surv_Sm_Ps_0.01_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.01_2)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30       30     15.8     12.77        24
## primary=S.m 30       26     40.2      5.02        24
## 
##  Chisq= 24  on 1 degrees of freedom, p= 9e-07
surv_Sm_Ps_0.1_2<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.1_2)
surv_Sm_Ps_0.1_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.1_2)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30       30     14.6     16.35      40.9
## primary=S.m 30       30     45.4      5.24      40.9
## 
##  Chisq= 40.9  on 1 degrees of freedom, p= 2e-10
surv_Sm_Ps_1.0_2<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_1.0_2)
surv_Sm_Ps_1.0_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_1.0_2)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30       30     13.5     20.14        44
## primary=S.m 30       30     46.5      5.85        44
## 
##  Chisq= 44  on 1 degrees of freedom, p= 3e-11
plot(survfit(Surv(Sm_Ps_block2$hour, Sm_Ps_block2$status) ~ Sm_Ps_block2$primary + Sm_Ps_block2$secondary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=2.3, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)

legend("bottomleft", bty="n" , c("PBS-PBS", "0.001", "PBS-0.01", "PBS-0.1", "PBS-1.0", "S.m-PBS", "S.m-0.001", "S.m-0.01", "S.m-0.1", "S.m-1.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=4, cex=1)

Block 3

Sm_Ps_block3 <- subset(Sm_Ps, date == "10/24/2020")
table(Sm_Ps_block3$primary, Sm_Ps_block3$secondary_dose)
##      
##        0 0.001 0.01 0.1  1
##   PBS 20    22   28  30 15
##   S.m 19    20   26  22 26
#subset by secondary dose
Sm_Ps_0.001_3<-subset(Sm_Ps_block3, secondary_dose=="0.001")
Sm_Ps_0.01_3<-subset(Sm_Ps_block3, secondary_dose=="0.01")
Sm_Ps_0.1_3<-subset(Sm_Ps_block3, secondary_dose=="0.1")
Sm_Ps_1.0_3<-subset(Sm_Ps_block3, secondary_dose=="1")

surv_Sm_Ps_0.001_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.001_3)
surv_Sm_Ps_0.001_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.001_3)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 22       21     15.6      1.90      9.75
## primary=S.m 20       13     18.4      1.61      9.75
## 
##  Chisq= 9.8  on 1 degrees of freedom, p= 0.002
surv_Sm_Ps_0.01_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.01_3)
surv_Sm_Ps_0.01_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.01_3)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28       28     13.1     16.80      34.7
## primary=S.m 26       21     35.9      6.16      34.7
## 
##  Chisq= 34.7  on 1 degrees of freedom, p= 4e-09
surv_Sm_Ps_0.1_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.1_3)
surv_Sm_Ps_0.1_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.1_3)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30       30       14     18.33      40.2
## primary=S.m 22       21       37      6.93      40.2
## 
##  Chisq= 40.2  on 1 degrees of freedom, p= 2e-10
surv_Sm_Ps_1.0_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_1.0_3)
surv_Sm_Ps_1.0_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_1.0_3)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 15       15      4.8     21.71      34.2
## primary=S.m 26       26     36.2      2.88      34.2
## 
##  Chisq= 34.2  on 1 degrees of freedom, p= 5e-09
plot(survfit(Surv(Sm_Ps_block3$hour, Sm_Ps_block3$status) ~ Sm_Ps_block3$primary + Sm_Ps_block3$secondary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=2.3, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)

legend("bottomleft", bty="n" , c("PBS-PBS", "0.001", "PBS-0.01", "PBS-0.1", "PBS-1.0", "S.m-PBS", "S.m-0.001", "S.m-0.01", "S.m-0.1", "S.m-1.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=4, cex=1)

Statistical Tests

  • In order test whether there was protection provided by chronic infection for each secondary infectious dose, p-values from the log rank test for each date was combined using the Fisher’s test (experiments on all dates showed higher survival in conditions with chronic infection). Final p-value is reported in the text.
# Grab p-values from each individual logrank test (results shown above in each block)
pvalue_Sm_Ps_0.001_1<-pchisq(surv_Sm_Ps_0.001_1$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_0.01_1<-pchisq(surv_Sm_Ps_0.01_1$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_0.1_1<-pchisq(surv_Sm_Ps_0.1_1$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_1.0_1<-pchisq(surv_Sm_Ps_1.0_1$chisq, df=1, lower.tail=F)

pvalue_Sm_Ps_0.001_2<-pchisq(surv_Sm_Ps_0.001_2$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_0.01_2<-pchisq(surv_Sm_Ps_0.01_2$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_0.1_2<-pchisq(surv_Sm_Ps_0.1_2$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_1.0_2<-pchisq(surv_Sm_Ps_1.0_2$chisq, df=1, lower.tail=F)

pvalue_Sm_Ps_0.001_3<-pchisq(surv_Sm_Ps_0.001_3$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_0.01_3<-pchisq(surv_Sm_Ps_0.01_3$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_0.1_3<-pchisq(surv_Sm_Ps_0.1_3$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_1.0_3<-pchisq(surv_Sm_Ps_1.0_3$chisq, df=1, lower.tail=F)
# Use all p-values from a given condition to get Fisher's combined test p-value
pvalues_Sm_Ps_0.001 <- c(pvalue_Sm_Ps_0.001_2,pvalue_Sm_Ps_0.001_3)
test.stat_Sm_Ps_0.001 <- -2*sum(log(pvalues_Sm_Ps_0.001))
pfinal_Sm_Ps_0.001<-pchisq(test.stat_Sm_Ps_0.001, df=2*length(pvalues_Sm_Ps_0.001), lower.tail=F)

pvalues_Sm_Ps_0.01 <- c(pvalue_Sm_Ps_0.01_1,pvalue_Sm_Ps_0.01_2,pvalue_Sm_Ps_0.01_3)
test.stat_Sm_Ps_0.01 <- -2*sum(log(pvalues_Sm_Ps_0.01))
pfinal_Sm_Ps_0.01<-pchisq(test.stat_Sm_Ps_0.01, df=2*length(pvalues_Sm_Ps_0.01), lower.tail=F)


pvalues_Sm_Ps_0.1 <- c(pvalue_Sm_Ps_0.1_1,pvalue_Sm_Ps_0.1_2,pvalue_Sm_Ps_0.1_3)
test.stat_Sm_Ps_0.1 <- -2*sum(log(pvalues_Sm_Ps_0.1))
pfinal_Sm_Ps_0.1<-pchisq(test.stat_Sm_Ps_0.1, df=2*length(pvalues_Sm_Ps_0.1), lower.tail=F)

pvalues_Sm_Ps_1.0 <- c(pvalue_Sm_Ps_1.0_1,pvalue_Sm_Ps_1.0_2,pvalue_Sm_Ps_1.0_3)
test.stat_Sm_Ps_1.0 <- -2*sum(log(pvalues_Sm_Ps_1.0))
pfinal_Sm_Ps_1.0<-pchisq(test.stat_Sm_Ps_1.0, df=2*length(pvalues_Sm_Ps_1.0), lower.tail=F)
pvalues_Sm_Ps_0.001_table <- c(pvalue_Sm_Ps_0.001_1, pvalue_Sm_Ps_0.001_2, pvalue_Sm_Ps_0.001_3)
p_Sm_Ps_0.001<-scales::pvalue(pvalues_Sm_Ps_0.001_table)

pvalues_Sm_Ps_0.01_table <- c(pvalue_Sm_Ps_0.01_1,pvalue_Sm_Ps_0.01_2,pvalue_Sm_Ps_0.01_3)
p_Sm_Ps_0.01<-scales::pvalue(pvalues_Sm_Ps_0.01_table)

pvalues_Sm_Ps_0.1_table <- c(pvalue_Sm_Ps_0.1_1,pvalue_Sm_Ps_0.1_2,pvalue_Sm_Ps_0.1_3)
p_Sm_Ps_0.1<-scales::pvalue(pvalues_Sm_Ps_0.1_table)

pvalues_Sm_Ps_1.0_table <- c(pvalue_Sm_Ps_1.0_1,pvalue_Sm_Ps_1.0_2,pvalue_Sm_Ps_1.0_3)
p_Sm_Ps_1.0<-scales::pvalue(pvalues_Sm_Ps_1.0_table)



pvalue_table_Sm_Ps <- data.frame("0.001" = p_Sm_Ps_0.001, "0.01" = p_Sm_Ps_0.01, "0.1" = p_Sm_Ps_0.1, "1.0" = p_Sm_Ps_1.0)
colnames(pvalue_table_Sm_Ps) <- c(0.001,0.01, 0.1, 1.0)

pfisher_Sm_Ps<- c(pfinal_Sm_Ps_0.001, pfinal_Sm_Ps_0.01, pfinal_Sm_Ps_0.1, pfinal_Sm_Ps_1.0)
pfisher_Sm_Ps<-scales::pvalue(pfisher_Sm_Ps)
pvalue_table_Sm_Ps[nrow(pvalue_table_Sm_Ps) + 1,]<-pfisher_Sm_Ps

row.names(pvalue_table_Sm_Ps) <- c("Block 1 - October 1st, 2021", "Block 2 - October 10th, 2021", "Block 3 - October 24th, 2021", " Final P-value ")

pvalue_table_Sm_Ps %>%
  kbl(caption = "<center><strong>P-values for log-rank tests on individual doses within each experiment<center><strong>") %>%
   kable_classic(full_width = F, html_font = "Cambria", position = "center")%>%
    add_header_above(c(" "=1, "Infectious Dose (Abs600nm)" = 4)) %>%
    pack_rows("Individual Blocks", 1, 3) %>%
    pack_rows("Fisher's Combined", 4, 4)
P-values for log-rank tests on individual doses within each experiment
Infectious Dose (Abs600nm)
0.001 0.01 0.1 1
Individual Blocks
Block 1 - October 1st, 2021 0.287 <0.001 <0.001 <0.001
Block 2 - October 10th, 2021 0.001 <0.001 <0.001 <0.001
Block 3 - October 24th, 2021 0.002 <0.001 <0.001 <0.001
Fisher’s Combined
Final P-value <0.001 <0.001 <0.001 <0.001
  • Conclusion: Chronic infection with S. marcescens offered significat protection against all doses of P. sneebia (p<0.001)

Q4 Does the dose of the primary infection impact the strength of protection during secondary infection?

  • The infectious dose impacts the number of bacteria present in the fly during chronic infection (Chambers et al. 2019), therefore we hypothesized that this might influence the strength of protection during secondary infection.
  • Flies were injected with one of 5 doses of S. marcescens during primary infection (0.001, 0.01, 0.1, 1.0 and 2.0) or PBS. Then they were given a challenge dose of either P. rettgeri or P. sneebia (2.0 Abs600nM) as a secondary infection.

Study design: + Initially analyzed using Cox Proportional Hazards mixed effects model with date of injection, primary injection dose as fixed effects and random effect terms accounting for experimental structure (e.g. date of injection, housing of groups of flies in vials). However, this model violated the assumptions of Cox Proportional Hazards. + Therefore we switched to using log-rank pairwise comparisons to determine whether dose used to induce chronic infection significantly impacted survival of secondary infection by comparing to the control flies given only a secondary infection. + LogRank comparisons were done on data from each date individually + If experiments on all dates showed higher survival in conditions with chronic infection, p-values were combined using Fisher’s test and only the final p-value is reported in the text. + If there was a variable effect of chronic infection, the number of individual experiments with significant effects is reported. + Conclusions + Moderate doses of S. marcescens (Abs600nm = 0.01, 0.1 & 1.0) all provided significant protection against secondary infection with both P. rettgeri ad P. sneebia (P<0.001) + The lowest dose of S. marcescens (Abs600nm = 0.01, 0.1 & 1.0) did not provide protection against either secondary infection. + The highest dose provided inconsistent protection, sometimes providing robust protection and other times dying as quickly as (or even faster than) control flies carrying only a secondary infection.

P. rettgeri secondary infection

  • Conclusion: Moderate doses of the primary infection (Abs600nm = 0.01, 0.1 & 1.0) all provided significant protection against infection (P<0.001).
  • The lowest dose (Abs600nm = 0.001) did not provide significant protection
  • The highest dose provided inconsistent protection with 8 out of the 9 blocks demonstrating better survival than controls and one with worse survival than controls. In Block 1, flies given a primary infection at the highest dose died more than controls with only the secondary infection.

Figure 5A

Pr_varied_primary = subset(Pr_varied_primary, secondary == "P.rett")
Pr_varied_primary$primary_dose <- as.factor(Pr_varied_primary$primary_dose)
table(Pr_varied_primary$date, Pr_varied_primary$primary_dose)
##             
##               0 0.001 0.01 0.1  1  2
##   10/18/2020 40    38   35  40 27 31
##   10/20/2020 36    37   35  39 35 30
##   10/29/2020 40    37   46  36 37 33
##   11/2/2020  37    31   30  40 30 27
##   11/28/2020 39    47   47  48 45 42
##   11/7/2020  40    38   40  35 34 36
##   12/4/2020  35    38   36  34 27  0
##   6/17/2021  37    47   47  44 42 36
##   6/18/2021  33    48   48  42 42 63
table(Pr_varied_primary$primary_dose)
## 
##     0 0.001  0.01   0.1     1     2 
##   337   361   364   358   319   298
#png(file="Figure5A.png",  width = 1200, height = 900, res = 130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(Pr_varied_primary$day, Pr_varied_primary$status) ~ Pr_varied_primary$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i', ylim=c(0, 1), yaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=4.5)
box(lwd=5)

#legend("bottomleft", bty="n" , c("secondary infection only", "chronic infection (0.001)", "chronic infection (0.01)", "chronic infection (0.1)", "chronic infection (1.0)", "chronic infection (2.0)"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1.5)
#dev.off()

Individual blocks

Block 1

Pr_varied_primary_block1 <- subset(Pr_varied_primary, date == "10/18/2020")
table(Pr_varied_primary_block1$primary_dose)
## 
##     0 0.001  0.01   0.1     1     2 
##    40    38    35    40    27    31
#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_1<-subset(Pr_varied_primary_block1, primary_dose == "0.001"|primary_dose=="0")
table(Pr_varied_primary_0.001_1$primary_dose)
## 
##     0 0.001  0.01   0.1     1     2 
##    40    38     0     0     0     0
Pr_varied_primary_0.01_1<-subset(Pr_varied_primary_block1, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_1<-subset(Pr_varied_primary_block1, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_1<-subset(Pr_varied_primary_block1, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_1<-subset(Pr_varied_primary_block1, primary_dose=="2"|primary_dose=="0")

surv_Pr_varied_primary_0.001_1<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_1)
surv_Pr_varied_primary_0.001_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_1)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40       32     37.4     0.792      3.24
## primary=S.m 38       35     29.6     1.003      3.24
## 
##  Chisq= 3.2  on 1 degrees of freedom, p= 0.07
surv_Pr_varied_primary_0.01_1<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_1)
surv_Pr_varied_primary_0.01_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_1)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40       32     28.8     0.358      1.11
## primary=S.m 35       24     27.2     0.379      1.11
## 
##  Chisq= 1.1  on 1 degrees of freedom, p= 0.3
surv_Pr_varied_primary_0.1_1<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_1)
surv_Pr_varied_primary_0.1_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_1)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40       32     18.6      9.60      22.5
## primary=S.m 40       12     25.4      7.05      22.5
## 
##  Chisq= 22.5  on 1 degrees of freedom, p= 2e-06
surv_Pr_varied_primary_1.0_1<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_1)
surv_Pr_varied_primary_1.0_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_1)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40       32     26.3      1.22      3.27
## primary=S.m 27       22     27.7      1.16      3.27
## 
##  Chisq= 3.3  on 1 degrees of freedom, p= 0.07
surv_Pr_varied_primary_2.0_1<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_1)
surv_Pr_varied_primary_2.0_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_1)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40       32     37.4     0.776      4.16
## primary=S.m 31       31     25.6     1.133      4.16
## 
##  Chisq= 4.2  on 1 degrees of freedom, p= 0.04
par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block1$day, Pr_varied_primary_block1$status) ~ Pr_varied_primary_block1$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)

legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)

Block 2

Pr_varied_primary_block2 <- subset(Pr_varied_primary, date == "10/20/2020")
table(Pr_varied_primary_block2$primary_dose)
## 
##     0 0.001  0.01   0.1     1     2 
##    36    37    35    39    35    30
#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_2<-subset(Pr_varied_primary_block2, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_2<-subset(Pr_varied_primary_block2, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_2<-subset(Pr_varied_primary_block2, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_2<-subset(Pr_varied_primary_block2, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_2<-subset(Pr_varied_primary_block2, primary_dose=="2"|primary_dose=="0")

surv_Pr_varied_primary_0.001_2<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_2)
surv_Pr_varied_primary_0.001_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_2)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 36       31     32.6    0.0783     0.403
## primary=S.m 37       35     33.4    0.0764     0.403
## 
##  Chisq= 0.4  on 1 degrees of freedom, p= 0.5
surv_Pr_varied_primary_0.01_2<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_2)
surv_Pr_varied_primary_0.01_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_2)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 36       31     27.7     0.388      1.57
## primary=S.m 35       28     31.3     0.344      1.57
## 
##  Chisq= 1.6  on 1 degrees of freedom, p= 0.2
surv_Pr_varied_primary_0.1_2<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_2)
surv_Pr_varied_primary_0.1_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_2)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 36       31     19.6      6.63      16.7
## primary=S.m 39       18     29.4      4.42      16.7
## 
##  Chisq= 16.7  on 1 degrees of freedom, p= 4e-05
surv_Pr_varied_primary_1.0_2<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_2)
surv_Pr_varied_primary_1.0_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_2)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 36       31     21.7      4.03      9.17
## primary=S.m 35       28     37.3      2.34      9.17
## 
##  Chisq= 9.2  on 1 degrees of freedom, p= 0.002
surv_Pr_varied_primary_2.0_2<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_2)
surv_Pr_varied_primary_2.0_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_2)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 36       31     25.5     1.192      3.23
## primary=S.m 30       27     32.5     0.935      3.23
## 
##  Chisq= 3.2  on 1 degrees of freedom, p= 0.07
par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block2$day, Pr_varied_primary_block2$status) ~ Pr_varied_primary_block2$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)

legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)

Block 3

Pr_varied_primary_block3 <- subset(Pr_varied_primary, date == "10/29/2020")
table(Pr_varied_primary_block3$primary_dose)
## 
##     0 0.001  0.01   0.1     1     2 
##    40    37    46    36    37    33
#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_3<-subset(Pr_varied_primary_block3, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_3<-subset(Pr_varied_primary_block3, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_3<-subset(Pr_varied_primary_block3, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_3<-subset(Pr_varied_primary_block3, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_3<-subset(Pr_varied_primary_block3, primary_dose=="2"|primary_dose=="0")

surv_Pr_varied_primary_0.001_3<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_3)
surv_Pr_varied_primary_0.001_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_3)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40       36     34.6    0.0559     0.235
## primary=S.m 37       31     32.4    0.0597     0.235
## 
##  Chisq= 0.2  on 1 degrees of freedom, p= 0.6
surv_Pr_varied_primary_0.01_3<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_3)
surv_Pr_varied_primary_0.01_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_3)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40       36     28.6      1.93      5.61
## primary=S.m 46       32     39.4      1.40      5.61
## 
##  Chisq= 5.6  on 1 degrees of freedom, p= 0.02
surv_Pr_varied_primary_0.1_3<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_3)
surv_Pr_varied_primary_0.1_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_3)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40       36     18.3      17.1      41.4
## primary=S.m 36        9     26.7      11.7      41.4
## 
##  Chisq= 41.4  on 1 degrees of freedom, p= 1e-10
surv_Pr_varied_primary_1.0_3<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_3)
surv_Pr_varied_primary_1.0_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_3)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40       36     19.6     13.83      30.3
## primary=S.m 37       19     35.4      7.63      30.3
## 
##  Chisq= 30.3  on 1 degrees of freedom, p= 4e-08
surv_Pr_varied_primary_2.0_3<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_3)
surv_Pr_varied_primary_2.0_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_3)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40       36     27.8      2.43      6.33
## primary=S.m 33       27     35.2      1.92      6.33
## 
##  Chisq= 6.3  on 1 degrees of freedom, p= 0.01
par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block3$day, Pr_varied_primary_block3$status) ~ Pr_varied_primary_block3$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)

legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)

Block 4

Pr_varied_primary_block4 <- subset(Pr_varied_primary, date == "11/2/2020")
table(Pr_varied_primary_block4$primary_dose)
## 
##     0 0.001  0.01   0.1     1     2 
##    37    31    30    40    30    27
#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_4<-subset(Pr_varied_primary_block4, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_4<-subset(Pr_varied_primary_block4, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_4<-subset(Pr_varied_primary_block4, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_4<-subset(Pr_varied_primary_block4, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_4<-subset(Pr_varied_primary_block4, primary_dose=="2"|primary_dose=="0")

surv_Pr_varied_primary_0.001_4<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_4)
surv_Pr_varied_primary_0.001_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_4)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37       33     33.2  0.000817   0.00537
## primary=S.m 31       28     27.8  0.000973   0.00537
## 
##  Chisq= 0  on 1 degrees of freedom, p= 0.9
surv_Pr_varied_primary_0.01_4<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_4)
surv_Pr_varied_primary_0.01_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_4)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37       33     26.5      1.57      6.91
## primary=S.m 30       20     26.5      1.58      6.91
## 
##  Chisq= 6.9  on 1 degrees of freedom, p= 0.009
surv_Pr_varied_primary_0.1_4<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_4)
surv_Pr_varied_primary_0.1_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_4)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37       33     21.2      6.58      20.5
## primary=S.m 40       19     30.8      4.52      20.5
## 
##  Chisq= 20.5  on 1 degrees of freedom, p= 6e-06
surv_Pr_varied_primary_1.0_4<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_4)
surv_Pr_varied_primary_1.0_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_4)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37       33     20.9      7.08      19.2
## primary=S.m 30       19     31.1      4.74      19.2
## 
##  Chisq= 19.2  on 1 degrees of freedom, p= 1e-05
surv_Pr_varied_primary_2.0_4<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_4)
surv_Pr_varied_primary_2.0_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_4)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37       33     30.6     0.181     0.686
## primary=S.m 27       27     29.4     0.189     0.686
## 
##  Chisq= 0.7  on 1 degrees of freedom, p= 0.4
par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block4$day, Pr_varied_primary_block4$status) ~ Pr_varied_primary_block4$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)

legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)

Block 5

Pr_varied_primary_block5 <- subset(Pr_varied_primary, date == "11/7/2020")
table(Pr_varied_primary_block5$primary_dose)
## 
##     0 0.001  0.01   0.1     1     2 
##    40    38    40    35    34    36
#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_5<-subset(Pr_varied_primary_block5, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_5<-subset(Pr_varied_primary_block5, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_5<-subset(Pr_varied_primary_block5, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_5<-subset(Pr_varied_primary_block5, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_5<-subset(Pr_varied_primary_block5, primary_dose=="2"|primary_dose=="0")

surv_Pr_varied_primary_0.001_5<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_5)
surv_Pr_varied_primary_0.001_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_5)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40       35     33.9    0.0383     0.156
## primary=S.m 38       32     33.1    0.0391     0.156
## 
##  Chisq= 0.2  on 1 degrees of freedom, p= 0.7
surv_Pr_varied_primary_0.01_5<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_5)
surv_Pr_varied_primary_0.01_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_5)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40       35       32     0.289     0.942
## primary=S.m 40       35       38     0.243     0.942
## 
##  Chisq= 0.9  on 1 degrees of freedom, p= 0.3
surv_Pr_varied_primary_0.1_5<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_5)
surv_Pr_varied_primary_0.1_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_5)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40       35     21.5      8.41      21.9
## primary=S.m 35       15     28.5      6.37      21.9
## 
##  Chisq= 21.9  on 1 degrees of freedom, p= 3e-06
surv_Pr_varied_primary_1.0_5<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_5)
surv_Pr_varied_primary_1.0_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_5)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40       35       20     11.28      27.4
## primary=S.m 34       16       31      7.27      27.4
## 
##  Chisq= 27.4  on 1 degrees of freedom, p= 2e-07
surv_Pr_varied_primary_2.0_5<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_5)
surv_Pr_varied_primary_2.0_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_5)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40       35     26.1      3.02       6.7
## primary=S.m 36       34     42.9      1.84       6.7
## 
##  Chisq= 6.7  on 1 degrees of freedom, p= 0.01
par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block5$day, Pr_varied_primary_block5$status) ~ Pr_varied_primary_block5$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)

legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)

Block 6

Pr_varied_primary_block6 <- subset(Pr_varied_primary, date == "11/28/2020")
table(Pr_varied_primary_block6$primary_dose)
## 
##     0 0.001  0.01   0.1     1     2 
##    39    47    47    48    45    42
#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_6<-subset(Pr_varied_primary_block6, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_6<-subset(Pr_varied_primary_block6, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_6<-subset(Pr_varied_primary_block6, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_6<-subset(Pr_varied_primary_block6, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_6<-subset(Pr_varied_primary_block6, primary_dose=="2"|primary_dose=="0")

surv_Pr_varied_primary_0.001_6<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_6)
surv_Pr_varied_primary_0.001_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_6)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39       33     38.1     0.679      3.44
## primary=S.m 47       44     38.9     0.664      3.44
## 
##  Chisq= 3.4  on 1 degrees of freedom, p= 0.06
surv_Pr_varied_primary_0.01_6<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_6)
surv_Pr_varied_primary_0.01_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_6)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39       33     32.7   0.00346    0.0119
## primary=S.m 47       39     39.3   0.00288    0.0119
## 
##  Chisq= 0  on 1 degrees of freedom, p= 0.9
surv_Pr_varied_primary_0.1_6<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_6)
surv_Pr_varied_primary_0.1_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_6)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39       33     19.6      9.14      19.7
## primary=S.m 48       23     36.4      4.93      19.7
## 
##  Chisq= 19.7  on 1 degrees of freedom, p= 9e-06
surv_Pr_varied_primary_1.0_6<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_6)
surv_Pr_varied_primary_1.0_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_6)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39       33     19.1     10.03      21.9
## primary=S.m 45       20     33.9      5.67      21.9
## 
##  Chisq= 21.9  on 1 degrees of freedom, p= 3e-06
surv_Pr_varied_primary_2.0_6<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_6)
surv_Pr_varied_primary_2.0_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_6)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39       33     25.8      2.03      4.56
## primary=S.m 42       37     44.2      1.18      4.56
## 
##  Chisq= 4.6  on 1 degrees of freedom, p= 0.03
par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block6$day, Pr_varied_primary_block6$status) ~ Pr_varied_primary_block6$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)

legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)

Block 7

Pr_varied_primary_block7 <- subset(Pr_varied_primary, date == "12/4/2020")
table(Pr_varied_primary_block7$primary_dose)
## 
##     0 0.001  0.01   0.1     1     2 
##    35    38    36    34    27     0
#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_7<-subset(Pr_varied_primary_block7, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_7<-subset(Pr_varied_primary_block7, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_7<-subset(Pr_varied_primary_block7, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_7<-subset(Pr_varied_primary_block7, primary_dose=="1"|primary_dose=="0")
#Pr_varied_primary_2.0_7<-subset(Pr_varied_primary_block7, primary_dose=="2"|primary_dose=="0")

surv_Pr_varied_primary_0.001_7<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_7)
surv_Pr_varied_primary_0.001_7
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_7)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35       31     31.1  0.000230   0.00162
## primary=S.m 38       33     32.9  0.000217   0.00162
## 
##  Chisq= 0  on 1 degrees of freedom, p= 1
surv_Pr_varied_primary_0.01_7<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_7)
surv_Pr_varied_primary_0.01_7
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_7)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35       31     24.8      1.53      7.31
## primary=S.m 36       21     27.2      1.40      7.31
## 
##  Chisq= 7.3  on 1 degrees of freedom, p= 0.007
surv_Pr_varied_primary_0.1_7<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_7)
surv_Pr_varied_primary_0.1_7
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_7)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35       31     16.1      13.8      40.8
## primary=S.m 34        6     20.9      10.6      40.8
## 
##  Chisq= 40.8  on 1 degrees of freedom, p= 2e-10
surv_Pr_varied_primary_1.0_7<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_7)
surv_Pr_varied_primary_1.0_7
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_7)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35       31     20.8      4.96      14.1
## primary=S.m 27       21     31.2      3.31      14.1
## 
##  Chisq= 14.1  on 1 degrees of freedom, p= 2e-04
#surv_Pr_varied_primary_2.0_7<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_7)
#surv_Pr_varied_primary_2.0_7

par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block7$day, Pr_varied_primary_block7$status) ~ Pr_varied_primary_block7$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)

legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)

Block 8

Pr_varied_primary_block8 <- subset(Pr_varied_primary, date == "6/17/2021")
table(Pr_varied_primary_block8$primary_dose)
## 
##     0 0.001  0.01   0.1     1     2 
##    37    47    47    44    42    36
#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_8<-subset(Pr_varied_primary_block8, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_8<-subset(Pr_varied_primary_block8, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_8<-subset(Pr_varied_primary_block8, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_8<-subset(Pr_varied_primary_block8, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_8<-subset(Pr_varied_primary_block8, primary_dose=="2"|primary_dose=="0")

surv_Pr_varied_primary_0.001_8<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_8)
surv_Pr_varied_primary_0.001_8
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_8)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37       37     36.2    0.0185     0.123
## primary=S.m 47       44     44.8    0.0150     0.123
## 
##  Chisq= 0.1  on 1 degrees of freedom, p= 0.7
surv_Pr_varied_primary_0.01_8<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_8)
surv_Pr_varied_primary_0.01_8
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_8)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37       37     32.2     0.700      3.45
## primary=S.m 47       42     46.8     0.483      3.45
## 
##  Chisq= 3.4  on 1 degrees of freedom, p= 0.06
surv_Pr_varied_primary_0.1_8<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_8)
surv_Pr_varied_primary_0.1_8
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_8)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37       37     17.5      21.9      55.4
## primary=S.m 44       14     33.5      11.4      55.4
## 
##  Chisq= 55.4  on 1 degrees of freedom, p= 1e-13
surv_Pr_varied_primary_1.0_8<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_8)
surv_Pr_varied_primary_1.0_8
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_8)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37       37       16      27.4      67.1
## primary=S.m 42       13       34      12.9      67.1
## 
##  Chisq= 67.1  on 1 degrees of freedom, p= 3e-16
surv_Pr_varied_primary_2.0_8<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_8)
surv_Pr_varied_primary_2.0_8
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_8)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37       37       24      7.02      27.6
## primary=S.m 36       35       48      3.51      27.6
## 
##  Chisq= 27.6  on 1 degrees of freedom, p= 2e-07
par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block8$day, Pr_varied_primary_block8$status) ~ Pr_varied_primary_block8$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)

legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)

Block 9

Pr_varied_primary_block9<- subset(Pr_varied_primary, date == "6/18/2021")
table(Pr_varied_primary_block9$primary_dose)
## 
##     0 0.001  0.01   0.1     1     2 
##    33    48    48    42    42    63
#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_9<-subset(Pr_varied_primary_block9, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_9<-subset(Pr_varied_primary_block9, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_9<-subset(Pr_varied_primary_block9, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_9<-subset(Pr_varied_primary_block9, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_9<-subset(Pr_varied_primary_block9, primary_dose=="2"|primary_dose=="0")

surv_Pr_varied_primary_0.001_9<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_9)
surv_Pr_varied_primary_0.001_9
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_9)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 33       33     32.2    0.0206      1.39
## primary=S.m 48       46     46.8    0.0142      1.39
## 
##  Chisq= 1.4  on 1 degrees of freedom, p= 0.2
surv_Pr_varied_primary_0.01_9<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_9)
surv_Pr_varied_primary_0.01_9
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_9)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 33       33     27.7     1.013      10.5
## primary=S.m 48       38     43.3     0.648      10.5
## 
##  Chisq= 10.5  on 1 degrees of freedom, p= 0.001
surv_Pr_varied_primary_0.1_9<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_9)
surv_Pr_varied_primary_0.1_9
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_9)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 33       33     19.8      8.80      38.8
## primary=S.m 42       18     31.2      5.58      38.8
## 
##  Chisq= 38.8  on 1 degrees of freedom, p= 5e-10
surv_Pr_varied_primary_1.0_9<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_9)
surv_Pr_varied_primary_1.0_9
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_9)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 33       33     17.2      14.6      53.7
## primary=S.m 42       13     28.8       8.7      53.7
## 
##  Chisq= 53.7  on 1 degrees of freedom, p= 2e-13
surv_Pr_varied_primary_2.0_9<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_9)
surv_Pr_varied_primary_2.0_9
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_9)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 33       33     19.2      9.82      35.5
## primary=S.m 63       50     63.8      2.97      35.5
## 
##  Chisq= 35.5  on 1 degrees of freedom, p= 2e-09
par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block9$day, Pr_varied_primary_block9$status) ~ Pr_varied_primary_block9$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)

legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)

Statistical Tests

  • In order test whether there was protection provided by chronic infection for each secondary infectious dose, p-values from the log rank test for each date was combined using the Fisher’s test (experiments on all dates showed higher survival in conditions with chronic infection). Final p-value is reported in the text.
# Grab p-values from each individual logrank test (results shown above in each block)
pvalue_Pr_varied_primary_0.001_1<-pchisq(surv_Pr_varied_primary_0.001_1$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_1<-pchisq(surv_Pr_varied_primary_0.01_1$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_1<-pchisq(surv_Pr_varied_primary_0.1_1$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_1<-pchisq(surv_Pr_varied_primary_1.0_1$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_1<-pchisq(surv_Pr_varied_primary_2.0_1$chisq, df=1, lower.tail=F)

pvalue_Pr_varied_primary_0.001_2<-pchisq(surv_Pr_varied_primary_0.001_2$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_2<-pchisq(surv_Pr_varied_primary_0.01_2$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_2<-pchisq(surv_Pr_varied_primary_0.1_2$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_2<-pchisq(surv_Pr_varied_primary_1.0_2$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_2<-pchisq(surv_Pr_varied_primary_2.0_2$chisq, df=1, lower.tail=F)

pvalue_Pr_varied_primary_0.001_3<-pchisq(surv_Pr_varied_primary_0.001_3$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_3<-pchisq(surv_Pr_varied_primary_0.01_3$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_3<-pchisq(surv_Pr_varied_primary_0.1_3$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_3<-pchisq(surv_Pr_varied_primary_1.0_3$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_3<-pchisq(surv_Pr_varied_primary_2.0_3$chisq, df=1, lower.tail=F)

pvalue_Pr_varied_primary_0.001_4<-pchisq(surv_Pr_varied_primary_0.001_4$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_4<-pchisq(surv_Pr_varied_primary_0.01_4$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_4<-pchisq(surv_Pr_varied_primary_0.1_4$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_4<-pchisq(surv_Pr_varied_primary_1.0_4$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_4<-pchisq(surv_Pr_varied_primary_2.0_4$chisq, df=1, lower.tail=F)

pvalue_Pr_varied_primary_0.001_5<-pchisq(surv_Pr_varied_primary_0.001_5$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_5<-pchisq(surv_Pr_varied_primary_0.01_5$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_5<-pchisq(surv_Pr_varied_primary_0.1_5$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_5<-pchisq(surv_Pr_varied_primary_1.0_5$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_5<-pchisq(surv_Pr_varied_primary_2.0_5$chisq, df=1, lower.tail=F)

pvalue_Pr_varied_primary_0.001_6<-pchisq(surv_Pr_varied_primary_0.001_6$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_6<-pchisq(surv_Pr_varied_primary_0.01_6$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_6<-pchisq(surv_Pr_varied_primary_0.1_6$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_6<-pchisq(surv_Pr_varied_primary_1.0_6$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_6<-pchisq(surv_Pr_varied_primary_2.0_6$chisq, df=1, lower.tail=F)

pvalue_Pr_varied_primary_0.001_7<-pchisq(surv_Pr_varied_primary_0.001_7$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_7<-pchisq(surv_Pr_varied_primary_0.01_7$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_7<-pchisq(surv_Pr_varied_primary_0.1_7$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_7<-pchisq(surv_Pr_varied_primary_1.0_7$chisq, df=1, lower.tail=F)
#pvalue_Pr_varied_primary_2.0_7<-pchisq(surv_Pr_varied_primary_2.0_7$chisq, df=1, lower.tail=F)

pvalue_Pr_varied_primary_0.001_8<-pchisq(surv_Pr_varied_primary_0.001_8$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_8<-pchisq(surv_Pr_varied_primary_0.01_8$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_8<-pchisq(surv_Pr_varied_primary_0.1_8$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_8<-pchisq(surv_Pr_varied_primary_1.0_8$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_8<-pchisq(surv_Pr_varied_primary_2.0_8$chisq, df=1, lower.tail=F)

pvalue_Pr_varied_primary_0.001_9<-pchisq(surv_Pr_varied_primary_0.001_9$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_9<-pchisq(surv_Pr_varied_primary_0.01_9$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_9<-pchisq(surv_Pr_varied_primary_0.1_9$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_9<-pchisq(surv_Pr_varied_primary_1.0_9$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_9<-pchisq(surv_Pr_varied_primary_2.0_9$chisq, df=1, lower.tail=F)
# Use all p-values from a given condition to get Fisher's combined test p-value
pvalues_Pr_varied_primary_0.001 <- c(pvalue_Pr_varied_primary_0.001_2 ,pvalue_Pr_varied_primary_0.001_3, pvalue_Pr_varied_primary_0.001_4, pvalue_Pr_varied_primary_0.001_5, pvalue_Pr_varied_primary_0.001_6, pvalue_Pr_varied_primary_0.001_7, pvalue_Pr_varied_primary_0.001_8, pvalue_Pr_varied_primary_0.001_9)
test.stat_Pr_varied_primary_0.001 <- -2*sum(log(pvalues_Pr_varied_primary_0.001))
pfinal_Pr_varied_primary_0.001<-pchisq(test.stat_Pr_varied_primary_0.001, df=2*length(pvalues_Pr_varied_primary_0.001), lower.tail=F)

pvalues_Pr_varied_primary_0.01 <- c(pvalue_Pr_varied_primary_0.01_1,pvalue_Pr_varied_primary_0.01_2,pvalue_Pr_varied_primary_0.01_3,pvalue_Pr_varied_primary_0.01_4,pvalue_Pr_varied_primary_0.01_5,pvalue_Pr_varied_primary_0.01_6,pvalue_Pr_varied_primary_0.01_7,pvalue_Pr_varied_primary_0.01_8,pvalue_Pr_varied_primary_0.01_9)
test.stat_Pr_varied_primary_0.01 <- -2*sum(log(pvalues_Pr_varied_primary_0.01))
pfinal_Pr_varied_primary_0.01<-pchisq(test.stat_Pr_varied_primary_0.01, df=2*length(pvalues_Pr_varied_primary_0.01), lower.tail=F)


pvalues_Pr_varied_primary_0.1 <- c(pvalue_Pr_varied_primary_0.1_1,pvalue_Pr_varied_primary_0.1_2,pvalue_Pr_varied_primary_0.1_3,pvalue_Pr_varied_primary_0.1_4,pvalue_Pr_varied_primary_0.1_5,pvalue_Pr_varied_primary_0.1_6,pvalue_Pr_varied_primary_0.1_7,pvalue_Pr_varied_primary_0.1_8,pvalue_Pr_varied_primary_0.1_9)
test.stat_Pr_varied_primary_0.1 <- -2*sum(log(pvalues_Pr_varied_primary_0.1))
pfinal_Pr_varied_primary_0.1<-pchisq(test.stat_Pr_varied_primary_0.1, df=2*length(pvalues_Pr_varied_primary_0.1), lower.tail=F)

pvalues_Pr_varied_primary_1.0 <- c(pvalue_Pr_varied_primary_1.0_1,pvalue_Pr_varied_primary_1.0_2,pvalue_Pr_varied_primary_1.0_3,pvalue_Pr_varied_primary_1.0_4,pvalue_Pr_varied_primary_1.0_5,pvalue_Pr_varied_primary_1.0_6,pvalue_Pr_varied_primary_1.0_7,pvalue_Pr_varied_primary_1.0_8,pvalue_Pr_varied_primary_1.0_9)
test.stat_Pr_varied_primary_1.0 <- -2*sum(log(pvalues_Pr_varied_primary_1.0))
pfinal_Pr_varied_primary_1.0<-pchisq(test.stat_Pr_varied_primary_1.0, df=2*length(pvalues_Pr_varied_primary_1.0), lower.tail=F)

pvalues_Pr_varied_primary_2.0 <- c(pvalue_Pr_varied_primary_2.0_1,pvalue_Pr_varied_primary_2.0_2,pvalue_Pr_varied_primary_2.0_3,pvalue_Pr_varied_primary_2.0_4,pvalue_Pr_varied_primary_2.0_5,pvalue_Pr_varied_primary_2.0_6,pvalue_Pr_varied_primary_2.0_8,pvalue_Pr_varied_primary_2.0_9)
test.stat_Pr_varied_primary_2.0 <- -2*sum(log(pvalues_Pr_varied_primary_2.0))
pfinal_Pr_varied_primary_2.0<-pchisq(test.stat_Pr_varied_primary_2.0, df=2*length(pvalues_Pr_varied_primary_2.0), lower.tail=F)
pvalues_Pr_varied_primary_0.001_table <- c(pvalue_Pr_varied_primary_0.001_1, pvalue_Pr_varied_primary_0.001_2, pvalue_Pr_varied_primary_0.001_3, pvalue_Pr_varied_primary_0.001_4, pvalue_Pr_varied_primary_0.001_5, pvalue_Pr_varied_primary_0.001_6, pvalue_Pr_varied_primary_0.001_7, pvalue_Pr_varied_primary_0.001_8, pvalue_Pr_varied_primary_0.001_9)
p_Pr_varied_primary_0.001<-scales::pvalue(pvalues_Pr_varied_primary_0.001_table)

pvalues_Pr_varied_primary_0.01_table <- c(pvalue_Pr_varied_primary_0.01_1,pvalue_Pr_varied_primary_0.01_2,pvalue_Pr_varied_primary_0.01_3,pvalue_Pr_varied_primary_0.01_4,pvalue_Pr_varied_primary_0.01_5,pvalue_Pr_varied_primary_0.01_6,pvalue_Pr_varied_primary_0.01_7,pvalue_Pr_varied_primary_0.01_8,pvalue_Pr_varied_primary_0.01_9)
p_Pr_varied_primary_0.01<-scales::pvalue(pvalues_Pr_varied_primary_0.01_table)

pvalues_Pr_varied_primary_0.1_table <- c(pvalue_Pr_varied_primary_0.1_1,pvalue_Pr_varied_primary_0.1_2, pvalue_Pr_varied_primary_0.1_3, pvalue_Pr_varied_primary_0.1_4, pvalue_Pr_varied_primary_0.1_5, pvalue_Pr_varied_primary_0.1_6, pvalue_Pr_varied_primary_0.1_7, pvalue_Pr_varied_primary_0.1_8, pvalue_Pr_varied_primary_0.1_9)
p_Pr_varied_primary_0.1<-scales::pvalue(pvalues_Pr_varied_primary_0.1_table)

pvalues_Pr_varied_primary_1.0_table <- c(pvalue_Pr_varied_primary_1.0_1,pvalue_Pr_varied_primary_1.0_2, pvalue_Pr_varied_primary_1.0_3, pvalue_Pr_varied_primary_1.0_4, pvalue_Pr_varied_primary_1.0_5, pvalue_Pr_varied_primary_1.0_6, pvalue_Pr_varied_primary_1.0_7, pvalue_Pr_varied_primary_1.0_8, pvalue_Pr_varied_primary_1.0_9)
p_Pr_varied_primary_1.0<-scales::pvalue(pvalues_Pr_varied_primary_1.0_table)

pvalues_Pr_varied_primary_2.0_table <- c(pvalue_Pr_varied_primary_2.0_1,pvalue_Pr_varied_primary_2.0_2, pvalue_Pr_varied_primary_2.0_3, pvalue_Pr_varied_primary_2.0_4, pvalue_Pr_varied_primary_2.0_5, pvalue_Pr_varied_primary_2.0_6, NA, pvalue_Pr_varied_primary_2.0_8, pvalue_Pr_varied_primary_2.0_9)
p_Pr_varied_primary_2.0<-scales::pvalue(pvalues_Pr_varied_primary_2.0_table)
p_Pr_varied_primary_2.0 <- p_Pr_varied_primary_2.0 %>% replace_na("--")


pvalue_table_Pr_varied_primary <- data.frame("0.001" = p_Pr_varied_primary_0.001, "0.01" = p_Pr_varied_primary_0.01, "0.1" = p_Pr_varied_primary_0.1, "1.0" = p_Pr_varied_primary_1.0, "2.0" = p_Pr_varied_primary_2.0)
colnames(pvalue_table_Pr_varied_primary) <- c(0.001,0.01, 0.1, 1.0,2.0)

pfisher_Pr_varied_primary<- c(pfinal_Pr_varied_primary_0.001, pfinal_Pr_varied_primary_0.01, pfinal_Pr_varied_primary_0.1, pfinal_Pr_varied_primary_1.0, NA)
pfisher_Pr_varied_primary<-scales::pvalue(pfisher_Pr_varied_primary)
pfisher_Pr_varied_primary<- pfisher_Pr_varied_primary%>% replace_na("N.C.")
pvalue_table_Pr_varied_primary[nrow(pvalue_table_Pr_varied_primary) + 1,]<-pfisher_Pr_varied_primary

row.names(pvalue_table_Pr_varied_primary) <- c("Block 1 - October 18th, 2020", "Block 2 - October 20th, 2020", "Block 3 - October 29th, 2020","Block 4 - November 2nd, 2020", "Block 5 - November 7th, 2020", "Block 6 - November 28th, 2020","Block 7 - December 4th, 2020", "Block 8 - June 17th, 2021", "Block 9 - June 18th, 2021", " Final P-value ")

pvalue_table_Pr_varied_primary %>%
  kbl(caption = "<center><strong>P-values for log-rank tests on individual doses within each experiment<center><strong>") %>%
   kable_classic(full_width = F, html_font = "Cambria", position = "center")%>%
    add_header_above(c(" "=1, "Primary Infectious Dose (Abs600nm)" = 5)) %>%
    pack_rows("Individual Blocks", 1, 9) %>%
    pack_rows("Fisher's Combined", 10, 10)
P-values for log-rank tests on individual doses within each experiment
Primary Infectious Dose (Abs600nm)
0.001 0.01 0.1 1 2
Individual Blocks
Block 1 - October 18th, 2020 0.072 0.293 <0.001 0.070 0.041
Block 2 - October 20th, 2020 0.526 0.211 <0.001 0.002 0.072
Block 3 - October 29th, 2020 0.628 0.018 <0.001 <0.001 0.012
Block 4 - November 2nd, 2020 0.942 0.009 <0.001 <0.001 0.407
Block 5 - November 7th, 2020 0.693 0.332 <0.001 <0.001 0.010
Block 6 - November 28th, 2020 0.064 0.913 <0.001 <0.001 0.033
Block 7 - December 4th, 2020 0.968 0.007 <0.001 <0.001 –
Block 8 - June 17th, 2021 0.726 0.063 <0.001 <0.001 <0.001
Block 9 - June 18th, 2021 0.238 0.001 <0.001 <0.001 <0.001
Fisher’s Combined
Final P-value 0.733 <0.001 <0.001 <0.001 N.C.
  • Conclusion: Moderate doses of the primary infection (Abs600nm = 0.01, 0.1 & 1.0) all provided significant protection against infection (P<0.001). The lowest dose (Abs600nm = 0.001) did not provide significant protection (p=0.76). The highest dose provided inconsistent protection with 8 out of the 9 replicates demonstrating better survival than controls and one with worse survival than controls therefore we could not calculate a combine p-value.

P. sneebia secondary infection

  • Conclusion: All doses used to induce primary infection provided significant protection (P < 0.001), except for the lowest dose (Abs600nm = 0.261). Note: the highest primary dose provided inconsistent protection (4 out of 5 blocks had significant protection for all others, for Block 1 the mortality was the same as flies carrying only secondary infection)

Figure 5B

Ps_varied_primary <- subset(Ps_varied_primary, secondary == "P.snee")
table(Ps_varied_primary$date, Ps_varied_primary$primary_dose)
##            
##              0 0.001 0.01 0.1  1  2
##   3/10/2021 38    49   48  49 30  6
##   3/17/2021 37    50   46  42 22  0
##   3/31/2021 39    46   47  42 43 52
##   4/21/2021 35    44   46  45 30 49
##   4/28/2021 37    45   41  43 38 34
##   6/20/2021 27    48   49  48 47 60
table(Ps_varied_primary$primary_dose)
## 
##     0 0.001  0.01   0.1     1     2 
##   213   282   277   269   210   201
#png(file="Figure5B.png",  width = 1200, height = 900, res = 130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(Ps_varied_primary$hour, Ps_varied_primary$status) ~ Ps_varied_primary$condition), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="hours post-secondary infection", cex.lab=2.5, xlim=c(0,48), xaxs='i', ylim=c(0, 1), yaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)

legend("bottomleft", bty="n" , c("secondary infection only", "chronic (0.001) & secondary infection", "chronic (0.01) & secondary infection", "chronic (0.1) & secondary infection", "chronic (1.0) & secondary infection", "chronic (2.0) & secondary infection"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1.5)

#dev.off()

Individual Blocks

Block 1

Ps_varied_primary_block_1 <- subset(Ps_varied_primary, date == "3/10/2021")
table(Ps_varied_primary_block_1$primary_dose)
## 
##     0 0.001  0.01   0.1     1     2 
##    38    49    48    49    30     6
#subset by primary dose for pairwise comparison dose
Ps_varied_primary_0.001_1<-subset(Ps_varied_primary_block_1, primary_dose == "0.001"|primary_dose=="0")
table(Ps_varied_primary_0.001_1$primary_dose)
## 
##     0 0.001 
##    38    49
Ps_varied_primary_0.01_1<-subset(Ps_varied_primary_block_1, primary_dose=="0.01"|primary_dose=="0")
Ps_varied_primary_0.1_1<-subset(Ps_varied_primary_block_1, primary_dose=="0.1"|primary_dose=="0")
Ps_varied_primary_1.0_1<-subset(Ps_varied_primary_block_1, primary_dose=="1"|primary_dose=="0")
Ps_varied_primary_2.0_1<-subset(Ps_varied_primary_block_1, primary_dose=="2"|primary_dose=="0")

surv_Ps_varied_primary_0.001_1<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.001_1)
surv_Ps_varied_primary_0.001_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.001_1)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 38       38     30.3      1.98      4.76
## primary=S.m 49       47     54.7      1.10      4.76
## 
##  Chisq= 4.8  on 1 degrees of freedom, p= 0.03
surv_Ps_varied_primary_0.01_1<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.01_1)
surv_Ps_varied_primary_0.01_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.01_1)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 38       38     31.6      1.29      3.14
## primary=S.m 48       48     54.4      0.75      3.14
## 
##  Chisq= 3.1  on 1 degrees of freedom, p= 0.08
surv_Ps_varied_primary_0.1_1<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.1_1)
surv_Ps_varied_primary_0.1_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.1_1)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 38       38     22.1     11.40      22.7
## primary=S.m 49       49     64.9      3.89      22.7
## 
##  Chisq= 22.7  on 1 degrees of freedom, p= 2e-06
surv_Ps_varied_primary_1.0_1<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_1.0_1)
surv_Ps_varied_primary_1.0_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_1.0_1)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 38       38     36.7    0.0435     0.148
## primary=S.m 30       30     31.3    0.0511     0.148
## 
##  Chisq= 0.1  on 1 degrees of freedom, p= 0.7
surv_Ps_varied_primary_2.0_1<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_2.0_1)
surv_Ps_varied_primary_2.0_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_2.0_1)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 38       38    37.12    0.0208     0.234
## primary=S.m  6        5     5.88    0.1314     0.234
## 
##  Chisq= 0.2  on 1 degrees of freedom, p= 0.6
plot(survfit(Surv(Ps_varied_primary_block_1$hour, Ps_varied_primary_block_1$status) ~ Ps_varied_primary_block_1$condition), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)

legend("bottomleft", bty="n" , c("PBS-P.snee 2.0", "S.m 0.001-P.snee 2.0", "S.m 0.01-P.snee 2.0", "S.m 0.1-P.snee 2.0", "S.m 1.0-P.snee 2.0", "S.m 2.0-P.snee 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)

Block 2

Ps_varied_primary_block_2 <- subset(Ps_varied_primary, date == "3/17/2021")
table(Ps_varied_primary_block_2$primary_dose)
## 
##     0 0.001  0.01   0.1     1 
##    37    50    46    42    22
#subset by primary dose for pairwise comparison dose
Ps_varied_primary_0.001_2<-subset(Ps_varied_primary_block_2, primary_dose == "0.001"|primary_dose=="0")
table(Ps_varied_primary_0.001_2$primary_dose)
## 
##     0 0.001 
##    37    50
Ps_varied_primary_0.01_2<-subset(Ps_varied_primary_block_2, primary_dose=="0.01"|primary_dose=="0")
Ps_varied_primary_0.1_2<-subset(Ps_varied_primary_block_2, primary_dose=="0.1"|primary_dose=="0")
Ps_varied_primary_1.0_2<-subset(Ps_varied_primary_block_2, primary_dose=="1"|primary_dose=="0")
#Ps_varied_primary_2.0_2<-subset(Ps_varied_primary_block_2, primary_dose=="2"|primary_dose=="0")

surv_Ps_varied_primary_0.001_2<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.001_2)
surv_Ps_varied_primary_0.001_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.001_2)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37       37     38.5    0.0568      0.16
## primary=S.m 50       49     47.5    0.0460      0.16
## 
##  Chisq= 0.2  on 1 degrees of freedom, p= 0.7
surv_Ps_varied_primary_0.01_2<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.01_2)
surv_Ps_varied_primary_0.01_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.01_2)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37       37       30     1.638      3.92
## primary=S.m 46       44       51     0.963      3.92
## 
##  Chisq= 3.9  on 1 degrees of freedom, p= 0.05
surv_Ps_varied_primary_0.1_2<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.1_2)
surv_Ps_varied_primary_0.1_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.1_2)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37       37     20.7     12.89      25.7
## primary=S.m 42       32     48.3      5.52      25.7
## 
##  Chisq= 25.7  on 1 degrees of freedom, p= 4e-07
surv_Ps_varied_primary_1.0_2<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_1.0_2)
surv_Ps_varied_primary_1.0_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_1.0_2)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37       37     27.8      3.04      9.36
## primary=S.m 22       14     23.2      3.64      9.36
## 
##  Chisq= 9.4  on 1 degrees of freedom, p= 0.002
#surv_Ps_varied_primary_2.0_2<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_2.0_2)
#surv_Ps_varied_primary_2.0_2

plot(survfit(Surv(Ps_varied_primary_block_2$hour, Ps_varied_primary_block_2$status) ~ Ps_varied_primary_block_2$condition), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)

legend("bottomleft", bty="n" , c("PBS-P.snee 2.0", "S.m 0.001-P.snee 2.0", "S.m 0.01-P.snee 2.0", "S.m 0.1-P.snee 2.0", "S.m 1.0-P.snee 2.0", "S.m 2.0-P.snee 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)

Block 3

Ps_varied_primary_block_3 <- subset(Ps_varied_primary, date == "3/31/2021")
table(Ps_varied_primary_block_3$primary_dose)
## 
##     0 0.001  0.01   0.1     1     2 
##    39    46    47    42    43    52
#subset by primary dose for pairwise comparison dose
Ps_varied_primary_0.001_3<-subset(Ps_varied_primary_block_3, primary_dose == "0.001"|primary_dose=="0")
table(Ps_varied_primary_0.001_3$primary_dose)
## 
##     0 0.001 
##    39    46
Ps_varied_primary_0.01_3<-subset(Ps_varied_primary_block_3, primary_dose=="0.01"|primary_dose=="0")
Ps_varied_primary_0.1_3<-subset(Ps_varied_primary_block_3, primary_dose=="0.1"|primary_dose=="0")
Ps_varied_primary_1.0_3<-subset(Ps_varied_primary_block_3, primary_dose=="1"|primary_dose=="0")
Ps_varied_primary_2.0_3<-subset(Ps_varied_primary_block_3, primary_dose=="2"|primary_dose=="0")

surv_Ps_varied_primary_0.001_3<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.001_3)
surv_Ps_varied_primary_0.001_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.001_3)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39       21       23     0.177     0.383
## primary=S.m 46       29       27     0.151     0.383
## 
##  Chisq= 0.4  on 1 degrees of freedom, p= 0.5
surv_Ps_varied_primary_0.01_3<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.01_3)
surv_Ps_varied_primary_0.01_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.01_3)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39       21     20.4    0.0196    0.0399
## primary=S.m 47       26     26.6    0.0150    0.0399
## 
##  Chisq= 0  on 1 degrees of freedom, p= 0.8
surv_Ps_varied_primary_0.1_3<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.1_3)
surv_Ps_varied_primary_0.1_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.1_3)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39       21     23.3     0.222     0.504
## primary=S.m 42       30     27.7     0.186     0.504
## 
##  Chisq= 0.5  on 1 degrees of freedom, p= 0.5
surv_Ps_varied_primary_1.0_3<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_1.0_3)
surv_Ps_varied_primary_1.0_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_1.0_3)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39       21     12.3      6.08      11.3
## primary=S.m 43        8     16.7      4.50      11.3
## 
##  Chisq= 11.3  on 1 degrees of freedom, p= 8e-04
surv_Ps_varied_primary_2.0_3<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_2.0_3)
surv_Ps_varied_primary_2.0_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_2.0_3)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39       21     10.7      9.96      17.2
## primary=S.m 52        7     17.3      6.15      17.2
## 
##  Chisq= 17.2  on 1 degrees of freedom, p= 3e-05
plot(survfit(Surv(Ps_varied_primary_block_3$hour, Ps_varied_primary_block_3$status) ~ Ps_varied_primary_block_3$condition), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)

legend("bottomleft", bty="n" , c("PBS-P.snee 2.0", "S.m 0.001-P.snee 2.0", "S.m 0.01-P.snee 2.0", "S.m 0.1-P.snee 2.0", "S.m 1.0-P.snee 2.0", "S.m 2.0-P.snee 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)

Block 4

Ps_varied_primary_block_4 <- subset(Ps_varied_primary, date == "4/21/2021")
table(Ps_varied_primary_block_4$primary_dose)
## 
##     0 0.001  0.01   0.1     1     2 
##    35    44    46    45    30    49
#subset by primary dose for pairwise comparison dose
Ps_varied_primary_0.001_4<-subset(Ps_varied_primary_block_4, primary_dose == "0.001"|primary_dose=="0")
table(Ps_varied_primary_0.001_4$primary_dose)
## 
##     0 0.001 
##    35    44
Ps_varied_primary_0.01_4<-subset(Ps_varied_primary_block_4, primary_dose=="0.01"|primary_dose=="0")
Ps_varied_primary_0.1_4<-subset(Ps_varied_primary_block_4, primary_dose=="0.1"|primary_dose=="0")
Ps_varied_primary_1.0_4<-subset(Ps_varied_primary_block_4, primary_dose=="1"|primary_dose=="0")
Ps_varied_primary_2.0_4<-subset(Ps_varied_primary_block_4, primary_dose=="2"|primary_dose=="0")

surv_Ps_varied_primary_0.001_4<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.001_4)
surv_Ps_varied_primary_0.001_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.001_4)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35       35     26.4      2.83      6.39
## primary=S.m 44       44     52.6      1.42      6.39
## 
##  Chisq= 6.4  on 1 degrees of freedom, p= 0.01
surv_Ps_varied_primary_0.01_4<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.01_4)
surv_Ps_varied_primary_0.01_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.01_4)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35       35     18.9     13.74        26
## primary=S.m 46       39     55.1      4.71        26
## 
##  Chisq= 26  on 1 degrees of freedom, p= 3e-07
surv_Ps_varied_primary_0.1_4<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.1_4)
surv_Ps_varied_primary_0.1_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.1_4)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35       35     12.4      41.3      69.8
## primary=S.m 45       22     44.6      11.5      69.8
## 
##  Chisq= 69.8  on 1 degrees of freedom, p= <2e-16
surv_Ps_varied_primary_1.0_4<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_1.0_4)
surv_Ps_varied_primary_1.0_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_1.0_4)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35       35       16      22.5      46.5
## primary=S.m 30        7       26      13.9      46.5
## 
##  Chisq= 46.5  on 1 degrees of freedom, p= 9e-12
surv_Ps_varied_primary_2.0_4<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_2.0_4)
surv_Ps_varied_primary_2.0_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_2.0_4)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35       35     19.5     12.33        23
## primary=S.m 49       26     41.5      5.79        23
## 
##  Chisq= 23  on 1 degrees of freedom, p= 2e-06
plot(survfit(Surv(Ps_varied_primary_block_4$hour, Ps_varied_primary_block_4$status) ~ Ps_varied_primary_block_4$condition), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)

legend("bottomleft", bty="n" , c("PBS-P.snee 2.0", "S.m 0.001-P.snee 2.0", "S.m 0.01-P.snee 2.0", "S.m 0.1-P.snee 2.0", "S.m 1.0-P.snee 2.0", "S.m 2.0-P.snee 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)

Block 5

Ps_varied_primary_block_5 <- subset(Ps_varied_primary, date == "4/28/2021")
table(Ps_varied_primary_block_5$primary_dose)
## 
##     0 0.001  0.01   0.1     1     2 
##    37    45    41    43    38    34
#subset by primary dose for pairwise comparison dose
Ps_varied_primary_0.001_5<-subset(Ps_varied_primary_block_5, primary_dose == "0.001"|primary_dose=="0")
table(Ps_varied_primary_0.001_5$primary_dose)
## 
##     0 0.001 
##    37    45
Ps_varied_primary_0.01_5<-subset(Ps_varied_primary_block_5, primary_dose=="0.01"|primary_dose=="0")
Ps_varied_primary_0.1_5<-subset(Ps_varied_primary_block_5, primary_dose=="0.1"|primary_dose=="0")
Ps_varied_primary_1.0_5<-subset(Ps_varied_primary_block_5, primary_dose=="1"|primary_dose=="0")
Ps_varied_primary_2.0_5<-subset(Ps_varied_primary_block_5, primary_dose=="2"|primary_dose=="0")

surv_Ps_varied_primary_0.001_5<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.001_5)
surv_Ps_varied_primary_0.001_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.001_5)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37       37     37.2  0.000842   0.00276
## primary=S.m 45       45     44.8  0.000699   0.00276
## 
##  Chisq= 0  on 1 degrees of freedom, p= 1
surv_Ps_varied_primary_0.01_5<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.01_5)
surv_Ps_varied_primary_0.01_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.01_5)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37       37     36.4   0.00838    0.0276
## primary=S.m 41       41     41.6   0.00735    0.0276
## 
##  Chisq= 0  on 1 degrees of freedom, p= 0.9
surv_Ps_varied_primary_0.1_5<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.1_5)
surv_Ps_varied_primary_0.1_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.1_5)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37       37     28.3      2.68       6.9
## primary=S.m 43       43     51.7      1.47       6.9
## 
##  Chisq= 6.9  on 1 degrees of freedom, p= 0.009
surv_Ps_varied_primary_1.0_5<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_1.0_5)
surv_Ps_varied_primary_1.0_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_1.0_5)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37       37     17.9     20.24      40.2
## primary=S.m 38       34     53.1      6.85      40.2
## 
##  Chisq= 40.2  on 1 degrees of freedom, p= 2e-10
surv_Ps_varied_primary_2.0_5<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_2.0_5)
surv_Ps_varied_primary_2.0_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_2.0_5)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37       37     18.8     17.71      37.4
## primary=S.m 34       30     48.2      6.89      37.4
## 
##  Chisq= 37.4  on 1 degrees of freedom, p= 9e-10
plot(survfit(Surv(Ps_varied_primary_block_5$hour, Ps_varied_primary_block_5$status) ~ Ps_varied_primary_block_5$condition), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)

legend("bottomleft", bty="n" , c("PBS-P.snee 2.0", "S.m 0.001-P.snee 2.0", "S.m 0.01-P.snee 2.0", "S.m 0.1-P.snee 2.0", "S.m 1.0-P.snee 2.0", "S.m 2.0-P.snee 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)

Block 6

Ps_varied_primary_block_6 <- subset(Ps_varied_primary, date == "6/20/2021")
table(Ps_varied_primary_block_6$primary_dose)
## 
##     0 0.001  0.01   0.1     1     2 
##    27    48    49    48    47    60
#subset by primary dose for pairwise comparison dose
Ps_varied_primary_0.001_6<-subset(Ps_varied_primary_block_6, primary_dose == "0.001"|primary_dose=="0")
table(Ps_varied_primary_0.001_6$primary_dose)
## 
##     0 0.001 
##    27    48
Ps_varied_primary_0.01_6<-subset(Ps_varied_primary_block_6, primary_dose=="0.01"|primary_dose=="0")
Ps_varied_primary_0.1_6<-subset(Ps_varied_primary_block_6, primary_dose=="0.1"|primary_dose=="0")
Ps_varied_primary_1.0_6<-subset(Ps_varied_primary_block_6, primary_dose=="1"|primary_dose=="0")
Ps_varied_primary_2.0_6<-subset(Ps_varied_primary_block_6, primary_dose=="2"|primary_dose=="0")

surv_Ps_varied_primary_0.001_6<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.001_6)
surv_Ps_varied_primary_0.001_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.001_6)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27       27     24.9    0.1730     0.436
## primary=S.m 48       48     50.1    0.0861     0.436
## 
##  Chisq= 0.4  on 1 degrees of freedom, p= 0.5
surv_Ps_varied_primary_0.01_6<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.01_6)
surv_Ps_varied_primary_0.01_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.01_6)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27       27     20.2     2.262      4.59
## primary=S.m 49       46     52.8     0.868      4.59
## 
##  Chisq= 4.6  on 1 degrees of freedom, p= 0.03
surv_Ps_varied_primary_0.1_6<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.1_6)
surv_Ps_varied_primary_0.1_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.1_6)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27       27     10.5     25.69      40.5
## primary=S.m 48       31     47.5      5.71      40.5
## 
##  Chisq= 40.5  on 1 degrees of freedom, p= 2e-10
surv_Ps_varied_primary_1.0_6<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_1.0_6)
surv_Ps_varied_primary_1.0_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_1.0_6)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27       27     7.82     47.01      71.3
## primary=S.m 47       19    38.18      9.63      71.3
## 
##  Chisq= 71.3  on 1 degrees of freedom, p= <2e-16
surv_Ps_varied_primary_2.0_6<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_2.0_6)
surv_Ps_varied_primary_2.0_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_2.0_6)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27       27     14.1     11.88        19
## primary=S.m 60       37     49.9      3.35        19
## 
##  Chisq= 19  on 1 degrees of freedom, p= 1e-05
plot(survfit(Surv(Ps_varied_primary_block_6$hour, Ps_varied_primary_block_6$status) ~ Ps_varied_primary_block_6$condition), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)

legend("bottomleft", bty="n" , c("PBS-P.snee 2.0", "S.m 0.001-P.snee 2.0", "S.m 0.01-P.snee 2.0", "S.m 0.1-P.snee 2.0", "S.m 1.0-P.snee 2.0", "S.m 2.0-P.snee 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)

Statistical Tests

  • In order test whether there was protection provided by chronic infection for each secondary infectious dose, p-values from the log rank test for each date was combined using the Fisher’s test (experiments on all dates showed higher survival in conditions with chronic infection). Final p-value is reported in the text.
# Grab p-values from each individual logrank test (results shown above in each block)
pvalue_Ps_varied_primary_0.001_1<-pchisq(surv_Ps_varied_primary_0.001_1$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.01_1<-pchisq(surv_Ps_varied_primary_0.01_1$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.1_1<-pchisq(surv_Ps_varied_primary_0.1_1$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_1.0_1<-pchisq(surv_Ps_varied_primary_1.0_1$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_2.0_1<-pchisq(surv_Ps_varied_primary_2.0_1$chisq, df=1, lower.tail=F)

pvalue_Ps_varied_primary_0.001_2<-pchisq(surv_Ps_varied_primary_0.001_2$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.01_2<-pchisq(surv_Ps_varied_primary_0.01_2$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.1_2<-pchisq(surv_Ps_varied_primary_0.1_2$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_1.0_2<-pchisq(surv_Ps_varied_primary_1.0_2$chisq, df=1, lower.tail=F)
#pvalue_Ps_varied_primary_2.0_2<-pchisq(surv_Ps_varied_primary_2.0_2$chisq, df=1, lower.tail=F)

pvalue_Ps_varied_primary_0.001_3<-pchisq(surv_Ps_varied_primary_0.001_3$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.01_3<-pchisq(surv_Ps_varied_primary_0.01_3$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.1_3<-pchisq(surv_Ps_varied_primary_0.1_3$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_1.0_3<-pchisq(surv_Ps_varied_primary_1.0_3$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_2.0_3<-pchisq(surv_Ps_varied_primary_2.0_3$chisq, df=1, lower.tail=F)

pvalue_Ps_varied_primary_0.001_4<-pchisq(surv_Ps_varied_primary_0.001_4$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.01_4<-pchisq(surv_Ps_varied_primary_0.01_4$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.1_4<-pchisq(surv_Ps_varied_primary_0.1_4$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_1.0_4<-pchisq(surv_Ps_varied_primary_1.0_4$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_2.0_4<-pchisq(surv_Ps_varied_primary_2.0_4$chisq, df=1, lower.tail=F)

pvalue_Ps_varied_primary_0.001_5<-pchisq(surv_Ps_varied_primary_0.001_5$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.01_5<-pchisq(surv_Ps_varied_primary_0.01_5$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.1_5<-pchisq(surv_Ps_varied_primary_0.1_5$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_1.0_5<-pchisq(surv_Ps_varied_primary_1.0_5$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_2.0_5<-pchisq(surv_Ps_varied_primary_2.0_5$chisq, df=1, lower.tail=F)

pvalue_Ps_varied_primary_0.001_6<-pchisq(surv_Ps_varied_primary_0.001_6$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.01_6<-pchisq(surv_Ps_varied_primary_0.01_6$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.1_6<-pchisq(surv_Ps_varied_primary_0.1_6$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_1.0_6<-pchisq(surv_Ps_varied_primary_1.0_6$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_2.0_6<-pchisq(surv_Ps_varied_primary_2.0_6$chisq, df=1, lower.tail=F)
# Use all p-values from a given condition to get Fisher's combined test p-value
pvalues_Ps_varied_primary_0.001 <- c(pvalue_Ps_varied_primary_0.001_2 ,pvalue_Ps_varied_primary_0.001_3, pvalue_Ps_varied_primary_0.001_4, pvalue_Ps_varied_primary_0.001_5, pvalue_Ps_varied_primary_0.001_6)
test.stat_Ps_varied_primary_0.001 <- -2*sum(log(pvalues_Ps_varied_primary_0.001))
pfinal_Ps_varied_primary_0.001<-pchisq(test.stat_Ps_varied_primary_0.001, df=2*length(pvalues_Ps_varied_primary_0.001), lower.tail=F)

pvalues_Ps_varied_primary_0.01 <- c(pvalue_Ps_varied_primary_0.01_1,pvalue_Ps_varied_primary_0.01_2,pvalue_Ps_varied_primary_0.01_3,pvalue_Ps_varied_primary_0.01_4,pvalue_Ps_varied_primary_0.01_5,pvalue_Ps_varied_primary_0.01_6)
test.stat_Ps_varied_primary_0.01 <- -2*sum(log(pvalues_Ps_varied_primary_0.01))
pfinal_Ps_varied_primary_0.01<-pchisq(test.stat_Ps_varied_primary_0.01, df=2*length(pvalues_Ps_varied_primary_0.01), lower.tail=F)


pvalues_Ps_varied_primary_0.1 <- c(pvalue_Ps_varied_primary_0.1_1,pvalue_Ps_varied_primary_0.1_2,pvalue_Ps_varied_primary_0.1_3,pvalue_Ps_varied_primary_0.1_4,pvalue_Ps_varied_primary_0.1_5,pvalue_Ps_varied_primary_0.1_6)
test.stat_Ps_varied_primary_0.1 <- -2*sum(log(pvalues_Ps_varied_primary_0.1))
pfinal_Ps_varied_primary_0.1<-pchisq(test.stat_Ps_varied_primary_0.1, df=2*length(pvalues_Ps_varied_primary_0.1), lower.tail=F)

pvalues_Ps_varied_primary_1.0 <- c(pvalue_Ps_varied_primary_1.0_1,pvalue_Ps_varied_primary_1.0_2,pvalue_Ps_varied_primary_1.0_3,pvalue_Ps_varied_primary_1.0_4,pvalue_Ps_varied_primary_1.0_5,pvalue_Ps_varied_primary_1.0_6)
test.stat_Ps_varied_primary_1.0 <- -2*sum(log(pvalues_Ps_varied_primary_1.0))
pfinal_Ps_varied_primary_1.0<-pchisq(test.stat_Ps_varied_primary_1.0, df=2*length(pvalues_Ps_varied_primary_1.0), lower.tail=F)

pvalues_Ps_varied_primary_2.0 <- c(pvalue_Ps_varied_primary_2.0_1,pvalue_Ps_varied_primary_2.0_3,pvalue_Ps_varied_primary_2.0_4,pvalue_Ps_varied_primary_2.0_5,pvalue_Ps_varied_primary_2.0_6)
test.stat_Ps_varied_primary_2.0 <- -2*sum(log(pvalues_Ps_varied_primary_2.0))
pfinal_Ps_varied_primary_2.0<-pchisq(test.stat_Ps_varied_primary_2.0, df=2*length(pvalues_Ps_varied_primary_2.0), lower.tail=F)
pvalues_Ps_varied_primary_0.001_table <- c(pvalue_Ps_varied_primary_0.001_1, pvalue_Ps_varied_primary_0.001_2, pvalue_Ps_varied_primary_0.001_3, pvalue_Ps_varied_primary_0.001_4, pvalue_Ps_varied_primary_0.001_5, pvalue_Ps_varied_primary_0.001_6)
p_Ps_varied_primary_0.001<-scales::pvalue(pvalues_Ps_varied_primary_0.001_table)

pvalues_Ps_varied_primary_0.01_table <- c(pvalue_Ps_varied_primary_0.01_1,pvalue_Ps_varied_primary_0.01_2,pvalue_Ps_varied_primary_0.01_3,pvalue_Ps_varied_primary_0.01_4,pvalue_Ps_varied_primary_0.01_5,pvalue_Ps_varied_primary_0.01_6)
p_Ps_varied_primary_0.01<-scales::pvalue(pvalues_Ps_varied_primary_0.01_table)

pvalues_Ps_varied_primary_0.1_table <- c(pvalue_Ps_varied_primary_0.1_1,pvalue_Ps_varied_primary_0.1_2, pvalue_Ps_varied_primary_0.1_3, pvalue_Ps_varied_primary_0.1_4, pvalue_Ps_varied_primary_0.1_5, pvalue_Ps_varied_primary_0.1_6)
p_Ps_varied_primary_0.1<-scales::pvalue(pvalues_Ps_varied_primary_0.1_table)

pvalues_Ps_varied_primary_1.0_table <- c(pvalue_Ps_varied_primary_1.0_1,pvalue_Ps_varied_primary_1.0_2, pvalue_Ps_varied_primary_1.0_3, pvalue_Ps_varied_primary_1.0_4, pvalue_Ps_varied_primary_1.0_5, pvalue_Ps_varied_primary_1.0_6)
p_Ps_varied_primary_1.0<-scales::pvalue(pvalues_Ps_varied_primary_1.0_table)

pvalues_Ps_varied_primary_2.0_table <- c(pvalue_Ps_varied_primary_2.0_1,NA, pvalue_Ps_varied_primary_2.0_3, pvalue_Ps_varied_primary_2.0_4, pvalue_Ps_varied_primary_2.0_5, pvalue_Ps_varied_primary_2.0_6)
p_Ps_varied_primary_2.0<-scales::pvalue(pvalues_Ps_varied_primary_2.0_table)
p_Ps_varied_primary_2.0 <- p_Ps_varied_primary_2.0 %>% replace_na("--")




pvalue_table_Ps_varied_primary <- data.frame("0.001" = p_Ps_varied_primary_0.001, "0.01" = p_Ps_varied_primary_0.01, "0.1" = p_Ps_varied_primary_0.1, "1.0" = p_Ps_varied_primary_1.0, "2.0" = p_Ps_varied_primary_2.0)
colnames(pvalue_table_Ps_varied_primary) <- c(0.001,0.01, 0.1, 1.0,2.0)

pfisher_Ps_varied_primary<- c(pfinal_Ps_varied_primary_0.001, pfinal_Ps_varied_primary_0.01, pfinal_Ps_varied_primary_0.1, pfinal_Ps_varied_primary_1.0, pfinal_Ps_varied_primary_2.0)
pfisher_Ps_varied_primary<-scales::pvalue(pfisher_Ps_varied_primary)
pvalue_table_Ps_varied_primary[nrow(pvalue_table_Ps_varied_primary) + 1,]<-pfisher_Ps_varied_primary

row.names(pvalue_table_Ps_varied_primary) <- c("Block 1 - March 10th, 2021", "Block 2 - March 17th, 2021", "Block 3 - March 31st, 2021","Block 4 - April 21t, 2021", "Block 5 - April 28th, 2021", "Block 6 - June 20th, 2021", " Final P-value ")

pvalue_table_Ps_varied_primary %>%
  kbl(caption = "<center><strong>P-values for log-rank tests on individual doses within each experiment<center><strong>") %>%
   kable_classic(full_width = F, html_font = "Cambria", position = "center")%>%
    add_header_above(c(" "=1, "Primary Infectious Dose (Abs600nm)" = 5)) %>%
    pack_rows("Individual Blocks", 1, 6) %>%
    pack_rows("Fisher's Combined", 7, 7)
P-values for log-rank tests on individual doses within each experiment
Primary Infectious Dose (Abs600nm)
0.001 0.01 0.1 1 2
Individual Blocks
Block 1 - March 10th, 2021 0.029 0.076 <0.001 0.700 0.628
Block 2 - March 17th, 2021 0.689 0.048 <0.001 0.002 –
Block 3 - March 31st, 2021 0.536 0.842 0.478 <0.001 <0.001
Block 4 - April 21t, 2021 0.011 <0.001 <0.001 <0.001 <0.001
Block 5 - April 28th, 2021 0.958 0.868 0.009 <0.001 <0.001
Block 6 - June 20th, 2021 0.509 0.032 <0.001 <0.001 <0.001
Fisher’s Combined
Final P-value 0.261 <0.001 <0.001 <0.001 <0.001
  • Conclusion: All doses used to induce primary infection provided significant protection (P < 0.001), except for the lowest dose (Abs600nm = 0.01)

Figure 5

q <- ggdraw() +
  draw_image("Figure5A.png")

r <- ggdraw() +
  draw_image("Figure5B.png")


multiplot5 <- q/r

v<-multiplot5 
 
  #plot_annotation(tag_levels = 'A')
 
ggsave("Fig.5.Primary.pdf", v, device = "pdf")
v